diff --git a/.github/workflows/test_build.yml b/.github/workflows/test_build.yml
new file mode 100644
index 0000000..0d3b2c7
--- /dev/null
+++ b/.github/workflows/test_build.yml
@@ -0,0 +1,103 @@
+name: mpecopt Test suite
+
+on:
+  pull_request:
+
+jobs:
+  # Main job that runs all the tests for mpecopt in matlab.
+  test_mpecopt:
+    strategy:
+      fail-fast: false # We always want to see where the failures are.
+      matrix:
+        matlab_version: ["R2021b", "R2022a", "R2022b", "R2023a", "R2023b", "R2024a"]
+        casadi:
+          - version: "3.5.5"
+            file: "casadi-linux-matlabR2014b-v3.5.5.tar.gz"
+            extract: "tar -zxvf"
+          - version: "3.6.7"
+            file: "casadi-3.6.7-linux64-matlab2018b.zip"
+            extract: "unzip"
+    runs-on: ubuntu-22.04
+
+    steps:
+    # Checkout repo
+    - name: Checkout mpecopt
+      uses: actions/checkout@v4
+      with:
+        submodules: 'recursive'
+    # Install matlab to run tests.
+    - name: Install MATLAB
+      uses: matlab-actions/setup-matlab@v2.3
+      with:
+        # NOTE: no longer earliest version to support more property validation functions
+        release: ${{ matrix.matlab_version }}
+        products: >
+          Parallel_Computing_Toolbox
+          Optimization_Toolbox
+        cache: true
+    # Download relevant matlab release artifacts for the casadi release we are using.
+    # TODO: Definitely cache this.
+    - name: Download CasADi release
+      uses: robinraju/release-downloader@v1.11
+      with:
+        repository: "casadi/casadi"
+        tag: ${{ matrix.casadi.version }}
+        fileName: ${{ matrix.casadi.file }}
+
+    # untar the casadi 3.5.5 release and then add it to the matlab path.
+    - name: Install CasADi
+      shell: bash
+      run: |
+        mkdir casadi
+        cd casadi
+        ${{ matrix.casadi.extract}} $GITHUB_WORKSPACE/${{ matrix.casadi.file }}
+        echo "MATLABPATH=$MATLABPATH:$GITHUB_WORKSPACE/casadi" >> $GITHUB_ENV
+
+    # Run matlab tests from the test directory and dump the results in junit format.
+    - name: Install mpecopt
+      uses: matlab-actions/run-command@v2.1
+      with:
+        command: addpath(genpath('src'));savepath;
+
+    # Run matlab tests from the test directory and dump the results in junit format.
+    - name: Run mpecopt tests
+      uses: matlab-actions/run-tests@v2.1
+      with:
+        source-folder: src
+        select-by-folder: test
+        test-results-junit: "test-results/matlab-results-${{ matrix.matlab_version }}-${{ matrix.casadi.version }}.xml"
+        code-coverage-cobertura: "test-results/coverage-${{ matrix.matlab_version }}-${{ matrix.casadi.version }}.xml"
+        strict: true # Non-explicitly suppressed warnings should fail.
+        use-parallel: true
+
+    # Run ReportGenerator to generage human readable coverage (only for latest versions of matlab and casadi for now)
+    - name: Setup .NET Core # Required to execute ReportGenerator
+      if: ${{ matrix.matlab_version == 'R2024a' && matrix.casadi.version == '3.6.7' }}
+      uses: actions/setup-dotnet@v4
+      with:
+        dotnet-version: 8.x
+        dotnet-quality: 'ga'
+    - name: ReportGenerator
+      if: ${{ matrix.matlab_version == 'R2024a' && matrix.casadi.version == '3.6.7' }}
+      uses: danielpalme/ReportGenerator-GitHub-Action@5.3.11
+      with:
+        reports: "test-results/coverage-${{ matrix.matlab_version }}-${{ matrix.casadi.version }}.xml"
+        targetdir: 'test-results'
+        reporttypes: 'HtmlInline;MarkdownSummaryGithub'
+        verbosity: 'Info' # The verbosity level of the log messages. Values: Verbose, Info, Warning, Error, Off
+        title: 'mpecopt coverage'
+    - name: Comment PR with coverage
+      if: ${{ matrix.matlab_version == 'R2024a' && matrix.casadi.version == '3.6.7' }}
+      uses: thollander/actions-comment-pull-request@v3
+      with:
+        github-token: ${{ secrets.GH_TOKEN }}
+        file-path: "test-results/SummaryGithub.md"
+        comment-tag: coverage
+        
+    # Upload raw coverage and test results.
+    - name: Upload test results
+      if: always()
+      uses: actions/upload-artifact@v4
+      with:
+        name: "test-results-${{ matrix.matlab_version }}-${{ matrix.casadi.version }}"
+        path: test-results/
\ No newline at end of file
diff --git a/src/Mpecopt.m b/src/Mpecopt.m
index d9ba40c..9271a9c 100644
--- a/src/Mpecopt.m
+++ b/src/Mpecopt.m
@@ -117,7 +117,9 @@
             h_comp_ii = full(mpec_casadi.h_comp_fun(x_k(1:dims.n_primal),solver_initialization.p0));
             h_std_ii = full(mpec_casadi.h_std_fun(x_k(1:dims.n_primal),solver_initialization.p0));
             f_opt_ii = full(mpec_casadi.f_fun(x_k(1:dims.n_primal),solver_initialization.p0));
-            print_iter_stats('I',0,f_opt_ii,h_std_ii,h_comp_ii,'/',0,'Initial guess',0,0,0,1)
+            if opts.verbose_solver
+                print_iter_stats('I',0,f_opt_ii,h_std_ii,h_comp_ii,'/',0,'Initial guess',0,0,0,1)
+            end
 
             t_phase_i_start = tic;
             switch opts.initialization_strategy
diff --git a/src/mpecopt_phase_ii.m b/src/mpecopt_phase_ii.m
index 7855055..251c294 100644
--- a/src/mpecopt_phase_ii.m
+++ b/src/mpecopt_phase_ii.m
@@ -361,7 +361,7 @@
 % --------------- compute multiplier-based stationary points --------------
 multiplier_based_stationarity_debug = stats.multiplier_based_stationarity;
 if (stats.success || k==settings.max_iter) && settings.compute_tnlp_stationary_point && phase_ii
-    if ~strcmp(settings.piece_nlp_strategy,'TNLP')
+    % if ~strcmp(settings.piece_nlp_strategy,'TNLP')
         % resolve TNLP for correct multipliers
         lbx_bnlp_k = lbx; ubx_bnlp_k = ubx;  % reset bounds of bnlp.
         lbg_tnlp_k = lbg; ubg_tnlp_k = ubg;
@@ -374,13 +374,13 @@
         ubx_bnlp_k(dims.ind_x1(active_set_estimate_k.I_00)) = 0;
         ubx_bnlp_k(dims.ind_x2(active_set_estimate_k.I_00)) = 0;
         settings.piece_nlp_strategy = initial_strategy;
+    % end
         t_nlp_start = tic;
         results_nlp = solver('x0',x_k,'p', p0, 'lbx', lbx_bnlp_k, 'ubx', ubx_bnlp_k,'lbg', lbg_tnlp_k, 'ubg', ubg_tnlp_k);
         cpu_time_nlp_k_l = toc(t_nlp_start);
         x_k_multi = full(results_nlp.x);
         lambda_x_k = full(results_nlp.lam_x);
-    end
-    [stats.multiplier_based_stationarity, ~] = determine_multipliers_based_stationary_point(x_k_multi,lambda_x_k,dims,settings);
+        [stats.multiplier_based_stationarity, ~] = determine_multipliers_based_stationary_point(x_k_multi,lambda_x_k,dims,settings);
 end
 
 % Debug falure of stationary point computation
diff --git a/test/TestPhaseOneOpts.m b/test/TestPhaseOneOpts.m
new file mode 100644
index 0000000..7324584
--- /dev/null
+++ b/test/TestPhaseOneOpts.m
@@ -0,0 +1,52 @@
+classdef TestPhaseOneOpts < matlab.unittest.TestCase
+    properties (TestParameter)
+        % Only test the 
+        initialization_strategy = {'RelaxAndProject', 'FeasibilityEll1General', 'FeasibilityEllInfGeneral', 'TakeInitialGuessActiveSet', 'TakeProvidedActiveSet'};
+    end
+
+    methods (Test, ParameterCombination = 'exhaustive')
+        function test_initialization_strategy(tc,initialization_strategy)
+            import casadi.*
+            import matlab.unittest.fixtures.SuppressedWarningsFixture
+            tc.applyFixture( ...
+                SuppressedWarningsFixture("optim:intlinprog:IgnoreX0"))
+
+            x1 = SX.sym('x1');
+            x2 = SX.sym('x2');
+            x3 = SX.sym('x3');
+            x4 = SX.sym('x4');
+            x5 = SX.sym('x5');
+            x6 = SX.sym('x6');
+            x7 = SX.sym('x7');
+            x8 = SX.sym('x8');
+
+            w = [x1;x2;x3;x4;x5;x6;x7;x8];
+
+            p = SX.sym('p');
+            f = (x1-5)^2+(2*x2+1)^2;
+            g = [2*(x2-1)-1.5*x1+x3-0.5*x4+x5
+                3*x1-x2-3-x6;...
+                -x1+0.5*x2+4-x7;...
+                -x1-x2+7-x8];
+            G = [x6;x7;x8];
+            H = [x3;x4;x5];
+            x0 = [0;0;1;0;0;0;1;1];
+            lbw = zeros(8,1);
+            ubw = inf*ones(8,1);
+
+            lbg = zeros(4,1);
+            ubg = zeros(4,1);
+
+            mpec = struct('x', w, 'f', f, 'g', g,'p',p,'G',G ,'H',H);
+            solver_initalization = struct('x0', x0, 'lbx',lbw, 'ubx',ubw,'lbg',lbg, 'ubg',ubg,'p0',1, 'y0', [0;1;1]);
+
+            opts = MPECOptimizerOptions();
+            opts.initialization_strategy = initialization_strategy;
+            %opts.verbose_solver = false;
+            %opts.verbose_summary = false;
+
+            solver = Mpecopt(mpec, opts);
+            [sol,stats] = solver.solve(solver_initalization);
+        end
+    end
+end
diff --git a/test/TestPhaseTwoOpts.m b/test/TestPhaseTwoOpts.m
new file mode 100644
index 0000000..60239aa
--- /dev/null
+++ b/test/TestPhaseTwoOpts.m
@@ -0,0 +1,52 @@
+classdef TestPhaseTwoOpts < matlab.unittest.TestCase
+    properties (TestParameter)
+        % Only test the 
+        piece_nlp_strategy = {'TNLP', 'BNLP_integer', 'BNLP_Gradient1', 'BNLP_Gradient2'};
+    end
+
+    methods (Test, ParameterCombination = 'exhaustive')
+        function test_initialization_strategy(tc,piece_nlp_strategy)
+            import casadi.*
+            import matlab.unittest.fixtures.SuppressedWarningsFixture
+            tc.applyFixture( ...
+                SuppressedWarningsFixture("optim:intlinprog:IgnoreX0"))
+
+            x1 = SX.sym('x1');
+            x2 = SX.sym('x2');
+            x3 = SX.sym('x3');
+            x4 = SX.sym('x4');
+            x5 = SX.sym('x5');
+            x6 = SX.sym('x6');
+            x7 = SX.sym('x7');
+            x8 = SX.sym('x8');
+
+            w = [x1;x2;x3;x4;x5;x6;x7;x8];
+
+            p = SX.sym('p');
+            f = (x1-5)^2+(2*x2+1)^2;
+            g = [2*(x2-1)-1.5*x1+x3-0.5*x4+x5
+                3*x1-x2-3-x6;...
+                -x1+0.5*x2+4-x7;...
+                -x1-x2+7-x8];
+            G = [x6;x7;x8];
+            H = [x3;x4;x5];
+            x0 = [0;0;1;0;0;0;1;1];
+            lbw = zeros(8,1);
+            ubw = inf*ones(8,1);
+
+            lbg = zeros(4,1);
+            ubg = zeros(4,1);
+
+            mpec = struct('x', w, 'f', f, 'g', g,'p',p,'G',G ,'H',H);
+            solver_initalization = struct('x0', x0, 'lbx',lbw, 'ubx',ubw,'lbg',lbg, 'ubg',ubg,'p0',1, 'y0', [0;1;1]);
+
+            opts = MPECOptimizerOptions();
+            opts.piece_nlp_strategy = piece_nlp_strategy;
+            %opts.verbose_solver = false;
+            %opts.verbose_summary = false;
+
+            solver = Mpecopt(mpec, opts);
+            [sol,stats] = solver.solve(solver_initalization);
+        end
+    end
+end