diff --git a/include/amici/defines.h b/include/amici/defines.h
index 3fcf43e58c..6f5dc6594f 100644
--- a/include/amici/defines.h
+++ b/include/amici/defines.h
@@ -87,8 +87,7 @@ typedef enum AMICI_sensi_order_TAG {
 typedef enum AMICI_sensi_meth_TAG {
     AMICI_SENSI_NONE,
     AMICI_SENSI_FSA,
-    AMICI_SENSI_ASA,
-    AMICI_SENSI_SS
+    AMICI_SENSI_ASA
 } AMICI_sensi_meth;
 
 /** linear solvers for CVODES/IDAS */
diff --git a/include/amici/solver.h b/include/amici/solver.h
index 80f997fb30..b7ed36e7bf 100644
--- a/include/amici/solver.h
+++ b/include/amici/solver.h
@@ -798,7 +798,7 @@ class Solver {
      * AMIFree frees allocation solver memory
      */
     virtual void AMIFree() = 0;
-
+    
     /**
      * AMIAdjInit initializes the adjoint problem
      *
@@ -807,7 +807,7 @@ class Solver {
      *
      */
     virtual void AMIAdjInit(long int steps, int interp) = 0;
-
+    
     /**
      * AMICreateB specifies solver method and initializes solver memory for the
      * backward problem
@@ -1122,7 +1122,6 @@ class Solver {
 
     /** flag indicating whether sensitivities are supposed to be computed */
     AMICI_sensi_order sensi = AMICI_SENSI_ORDER_NONE;
-
 };
 
 bool operator ==(const Solver &a, const Solver &b);
diff --git a/include/amici/solver_cvodes.h b/include/amici/solver_cvodes.h
index fee65ccae5..eaa424aabe 100644
--- a/include/amici/solver_cvodes.h
+++ b/include/amici/solver_cvodes.h
@@ -33,8 +33,7 @@ class CVodeSolver : public Solver {
      * @return The clone
      */
     virtual Solver* clone() const override;
-
-
+    
     void *AMICreate(int lmm, int iter) override;
 
     void AMISStolerances(double rtol, double atol) override;
@@ -95,6 +94,17 @@ class CVodeSolver : public Solver {
 
     int AMISolveF(realtype tout, AmiVector *yret, AmiVector *ypret, realtype *tret,
                   int itask, int *ncheckPtr) override;
+    
+
+    static int fxdot(realtype t, N_Vector x, N_Vector xdot, void *user_data);
+    
+    static int fJSparse(realtype t, N_Vector x, N_Vector xdot, SlsMat J,
+                        void *user_data, N_Vector tmp1, N_Vector tmp2,
+                        N_Vector tmp3);
+    
+    static int fJ(long int N, realtype t, N_Vector x, N_Vector xdot,
+                  DlsMat J, void *user_data, N_Vector tmp1,
+                  N_Vector tmp2, N_Vector tmp3);
 
     void AMISolveB(realtype tBout, int itaskB) override;
 
@@ -196,18 +206,10 @@ class CVodeSolver : public Solver {
 
     void setJacTimesVecFnB(int which) override;
     
-    static int fJ(long int N, realtype t, N_Vector x, N_Vector xdot,
-                    DlsMat J, void *user_data, N_Vector tmp1,
-                    N_Vector tmp2, N_Vector tmp3);
-    
     static int fJB(long int NeqBdot, realtype t, N_Vector x, N_Vector xB,
                    N_Vector xBdot, DlsMat JB, void *user_data, N_Vector tmp1B,
                    N_Vector tmp2B, N_Vector tmp3B);
     
-    static int fJSparse(realtype t, N_Vector x, N_Vector xdot, SlsMat J,
-                        void *user_data, N_Vector tmp1, N_Vector tmp2,
-                        N_Vector tmp3);
-    
     static int fJSparseB(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot,
                          SlsMat JB, void *user_data, N_Vector tmp1B,
                          N_Vector tmp2B, N_Vector tmp3B);
@@ -232,8 +234,6 @@ class CVodeSolver : public Solver {
     static int froot(realtype t, N_Vector x, realtype *root,
                      void *user_data);
     
-    static int fxdot(realtype t, N_Vector x, N_Vector xdot, void *user_data);
-    
     static int fxBdot(realtype t, N_Vector x, N_Vector xB,
                       N_Vector xBdot, void *user_data);
     
diff --git a/include/amici/steadystateproblem.h b/include/amici/steadystateproblem.h
index efca7bd27f..477ce25030 100644
--- a/include/amici/steadystateproblem.h
+++ b/include/amici/steadystateproblem.h
@@ -7,6 +7,8 @@
 
 #include <nvector/nvector_serial.h>
 
+#include <functional>
+#include <memory>
 
 namespace amici {
 
@@ -35,10 +37,11 @@ class SteadystateProblem {
                                 Model *model, int newton_status,
                                 double run_time, int it);
 
-    void getNewtonSimulation(ReturnData *rdata, Solver *solver,
-                                   Model *model, int it);
-    
+    void getSteadystateSimulation(ReturnData *rdata, Solver *solver,
+                                  Model *model, int it);
     
+    std::unique_ptr<void, std::function<void(void *)>> createSteadystateSimSolver(Solver *solver, Model *model, realtype tstart);
+
     /** default constructor
      * @param t pointer to time variable
      * @param x pointer to state variables
diff --git a/matlab/@amimodel/generateMatlabWrapper.m b/matlab/@amimodel/generateMatlabWrapper.m
index b59919db05..2778759509 100644
--- a/matlab/@amimodel/generateMatlabWrapper.m
+++ b/matlab/@amimodel/generateMatlabWrapper.m
@@ -312,19 +312,19 @@ function generateMatlabWrapper(nx, ny, np, nk, nz, o2flag, amimodelo2, wrapperFi
     fprintf(fid,'init = struct();\n');
     fprintf(fid,'if(~isempty(options_ami.x0))\n');
     fprintf(fid,'    if(size(options_ami.x0,2)~=1)\n');
-    fprintf(fid,'        error(''x0 field must be a row vector!'');\n');
+    fprintf(fid,'        error(''x0 field must be a column vector!'');\n');
     fprintf(fid,'    end\n');
     fprintf(fid,'    if(size(options_ami.x0,1)~=nxfull)\n');
-    fprintf(fid,'        error(''Number of columns in x0 field does not agree with number of states!'');\n');
+    fprintf(fid,'        error(''Number of rows in x0 field does not agree with number of states!'');\n');
     fprintf(fid,'    end\n');
     fprintf(fid,'    init.x0 = options_ami.x0;\n');
     fprintf(fid,'end\n');
     fprintf(fid,'if(~isempty(options_ami.sx0))\n');
     fprintf(fid,'    if(size(options_ami.sx0,2)~=nplist)\n');
-    fprintf(fid,'        error(''Number of rows in sx0 field does not agree with number of model parameters!'');\n');
+    fprintf(fid,'        error(''Number of columns in sx0 field does not agree with number of model parameters!'');\n');
     fprintf(fid,'    end\n');
     fprintf(fid,'    if(size(options_ami.sx0,1)~=nxfull)\n');
-    fprintf(fid,'        error(''Number of columns in sx0 field does not agree with number of states!'');\n');
+    fprintf(fid,'        error(''Number of rows in sx0 field does not agree with number of states!'');\n');
     fprintf(fid,'    end\n');
     fprintf(fid,'    init.sx0 = bsxfun(@times,options_ami.sx0,1./permute(chainRuleFactor(:),[2,1]));\n');
     fprintf(fid,'end\n');
diff --git a/python/amici/__init__.py b/python/amici/__init__.py
index d33cb3a0af..65b93200c5 100644
--- a/python/amici/__init__.py
+++ b/python/amici/__init__.py
@@ -120,7 +120,7 @@ def getFieldAsNumPyArray(rdata, field):
                        'x':   [rdata.nt, rdata.nx],
                        'x0':  [rdata.nx],
                        'sx':  [rdata.nt, rdata.nplist, rdata.nx],
-                       'sx0': [rdata.nx, rdata.nplist],
+                       'sx0': [rdata.nplist, rdata.nx],
                        
                        # observables
                        'y':       [rdata.nt, rdata.ny],
diff --git a/src/hdf5.cpp b/src/hdf5.cpp
index 52fcfaafa5..e06b7c642b 100644
--- a/src/hdf5.cpp
+++ b/src/hdf5.cpp
@@ -119,6 +119,14 @@ std::unique_ptr<ExpData> readSimulationExpData(std::string const& hdf5Filename,
         checkEventDimensionsCompatible(m, n, model);
     }
 
+    if(locationExists(file,  hdf5Root + "/condition")) {
+        edata->fixedParameters = getDoubleDataset1D(file, hdf5Root + "/condition");
+    }
+
+    if(locationExists(file,  hdf5Root + "/conditionPreequilibration")) {
+        edata->fixedParametersPreequilibration = getDoubleDataset1D(file, hdf5Root + "/conditionPreequilibration");
+    }
+
     return edata;
 }
 
@@ -171,7 +179,7 @@ void writeReturnData(const ReturnData &rdata, H5::H5File const& file, const std:
                                         rdata.nJ - 1, rdata.nplist);
 
     if (rdata.sx0.size())
-        createAndWriteDouble2DDataset(file, hdf5Location + "/sx0", rdata.sx0.data(), rdata.nx, rdata.nplist);
+        createAndWriteDouble2DDataset(file, hdf5Location + "/sx0", rdata.sx0.data(), rdata.nplist, rdata.nx);
 
     if (rdata.sx.size())
         createAndWriteDouble3DDataset(file, hdf5Location + "/sx", rdata.sx.data(), rdata.nt, rdata.nplist, rdata.nx);
@@ -499,7 +507,7 @@ void readModelDataFromHDF5(const H5::H5File &file, Model &model, const std::stri
         hsize_t length1 = 0;
         auto sx0 = getDoubleDataset2D(file, datasetPath + "/sx0", length0, length1);
         if(sx0.size()) {
-            if (length0 != (unsigned) model.nx && length1 != (unsigned) model.nplist())
+            if (length0 != (unsigned) model.nplist() && length1 != (unsigned) model.nx)
                 throw(AmiException("Dimension mismatch when reading sx0. Expected %dx%d, got %llu, %llu.",
                                    model.nx, model.nplist(), length0, length1));
             model.setInitialStateSensitivities(sx0);
diff --git a/src/returndata_matlab.cpp b/src/returndata_matlab.cpp
index a8c1279924..840ba9ebe2 100644
--- a/src/returndata_matlab.cpp
+++ b/src/returndata_matlab.cpp
@@ -52,6 +52,7 @@ mxArray *initMatlabReturnFields(ReturnData const *rdata) {
     mxArray *matlabSolutionStruct =
         mxCreateStructMatrix(1, 1, numFields, field_names_sol);
     
+    std::vector<int> perm0 = {1, 0};
     std::vector<int> perm1 = {0, 1};
     std::vector<int> perm2 = {0, 2, 1};
     std::vector<int> perm3 = {0, 2, 3, 1};
@@ -69,7 +70,7 @@ mxArray *initMatlabReturnFields(ReturnData const *rdata) {
     }
     if (rdata->nx > 0) {
         writeMatlabField2(matlabSolutionStruct, "x", rdata->x, rdata->nt, rdata->nx, perm1);
-        writeMatlabField2(matlabSolutionStruct, "x0",  rdata->x0, 1, rdata->nx, perm1);
+        writeMatlabField2(matlabSolutionStruct, "x0",  rdata->x0, rdata->nx, 1, perm1);
     }
     if (rdata->ny > 0) {
         writeMatlabField2(matlabSolutionStruct, "y", rdata->y, rdata->nt, rdata->ny, perm1);
@@ -77,7 +78,7 @@ mxArray *initMatlabReturnFields(ReturnData const *rdata) {
     }
     if (rdata->sensi >= AMICI_SENSI_ORDER_FIRST) {
         writeMatlabField1(matlabSolutionStruct, "sllh", rdata->sllh, rdata->nplist);
-        writeMatlabField2(matlabSolutionStruct, "sx0",  rdata->sx0, rdata->nx, rdata->nplist, perm1);
+        writeMatlabField2(matlabSolutionStruct, "sx0",  rdata->sx0, rdata->nplist, rdata->nx, perm0);
 
         if (rdata->sensi_meth == AMICI_SENSI_FSA) {
             writeMatlabField3(matlabSolutionStruct, "sx", rdata->sx, rdata->nt, rdata->nplist, rdata->nx, perm2);
diff --git a/src/solver.cpp b/src/solver.cpp
index 9a574bff4f..c5bf3555c5 100644
--- a/src/solver.cpp
+++ b/src/solver.cpp
@@ -86,12 +86,10 @@ void Solver::setupAMI(ForwardProblem *fwd, Model *model) {
             }
         }
 
-        if (sensi_meth == AMICI_SENSI_ASA) {
-            if (model->nx > 0) {
+        if (sensi_meth == AMICI_SENSI_ASA)
+            if (model->nx > 0)
                 /* Allocate space for the adjoint computation */
                 AMIAdjInit(maxsteps, interpType);
-            }
-        }
     }
 
     AMISetId(model);
diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp
index 1035da5a2a..55c822b36f 100644
--- a/src/steadystateproblem.cpp
+++ b/src/steadystateproblem.cpp
@@ -1,6 +1,7 @@
 #include "amici/steadystateproblem.h"
 #include "amici/model.h"
 #include "amici/solver.h"
+#include "amici/solver_cvodes.h"
 #include "amici/edata.h"
 #include "amici/forwardproblem.h"
 #include "amici/newton_solver.h"
@@ -11,6 +12,8 @@
 #include <ctime>
 #include <sundials/sundials_dense.h>
 #include <memory>
+#include <cvodes/cvodes.h>
+#include <cvodes/cvodes_klu.h>
 
 namespace amici {
 
@@ -48,7 +51,7 @@ void SteadystateProblem::workSteadyStateProblem(ReturnData *rdata,
     } catch(NewtonFailure& ex) {
         try {
             /* Newton solver did not find a steady state, so try integration */
-            getNewtonSimulation(rdata, solver, model, it);
+            getSteadystateSimulation(rdata, solver, model, it);
             newton_status = 2;
         } catch(AmiException& ex) {// may be integration failure from AmiSolve, so NewtonFailure won't do for all cases
             try {
@@ -57,29 +60,27 @@ void SteadystateProblem::workSteadyStateProblem(ReturnData *rdata,
             } catch(NewtonFailure& ex) {
                 // TODO: more informative NewtonFailure to give more informative error code
                 throw amici::IntegrationFailure(AMICI_CONV_FAILURE,*t);
-            } catch(...) {
-                throw AmiException("Internal error in steady state problem");
             }
-        } catch(...) {
-            throw AmiException("Internal error in steady state problem");
         }
     } catch(...) {
         throw AmiException("Internal error in steady state problem");
     }
     run_time = (double)((clock() - starttime) * 1000) / CLOCKS_PER_SEC;
-    getNewtonOutput(rdata, model, newton_status, run_time, it);
 
     /* Compute steady state sensitvities */
-    if (rdata->sensi_meth == AMICI_SENSI_FSA &&
-        rdata->sensi >= AMICI_SENSI_ORDER_FIRST)
+    if (rdata->sensi >= AMICI_SENSI_ORDER_FIRST && rdata->sensi_meth != AMICI_SENSI_NONE)
         newtonSolver.get()->getSensis(it, sx);
 
+    /* Get output of steady state solver, write it to x0 and reset time if necessary */
+    getNewtonOutput(rdata, model, newton_status, run_time, it);
+    
     /* Reinitialize solver with preequilibrated state */
     if (it == AMICI_PREEQUILIBRATE) {
         solver->AMIReInit(*t, x, &dx);
-        if (rdata->sensi >= AMICI_SENSI_ORDER_FIRST)
+        if (rdata->sensi >= AMICI_SENSI_ORDER_FIRST) {
             if (rdata->sensi_meth == AMICI_SENSI_FSA)
                 solver->AMISensReInit(solver->getInternalSensitivityMethod(), sx, &sdx);
+        }
     }
 }
 
@@ -231,8 +232,8 @@ void SteadystateProblem::getNewtonOutput(ReturnData *rdata,
 /* ----------------------------------------------------------------------------------
  */
 
-void SteadystateProblem::getNewtonSimulation(ReturnData *rdata, Solver *solver,
-                                            Model *model, int it) {
+void SteadystateProblem::getSteadystateSimulation(ReturnData *rdata, Solver *solver,
+                                                  Model *model, int it) {
     /**
      * Forward simulation is launched, if Newton solver fails in first try
      *
@@ -242,17 +243,21 @@ void SteadystateProblem::getNewtonSimulation(ReturnData *rdata, Solver *solver,
      * @param[in] it current timepoint index, <0 indicates preequilibration @type int
      */
  
-    realtype tstart;
-    
-    /* Newton solver did not work, so try a simulation: reinitialize solver */
-    if (it<1)
-        tstart = model->t0();
-    else
-        tstart = rdata->ts[it-1];
-    *t = tstart;
-    
-    model->fx0(x);
-    solver->AMIReInit(*t, x, &dx);
+    std::unique_ptr<void, std::function<void(void *)>> newton_sim;
+
+    /* Newton solver did not work, so try a simulation */
+    if (it<1) {
+        /* Preequilibration: Create a new CVode object for simulation */
+        *t = model->t0();
+        newton_sim = createSteadystateSimSolver(solver, model, *t);
+        model->fx0(x);
+    } else {
+        /* Carry on simulating from last point */
+        *t = rdata->ts[it-1];
+        model->fx0(x);
+        /* Reinitialize old solver */
+        solver->AMIReInit(*t, x, &dx);
+    }
     
     /* Loop over steps and check for convergence */
     double res_abs = INFINITY;
@@ -266,8 +271,14 @@ void SteadystateProblem::getNewtonSimulation(ReturnData *rdata, Solver *solver,
          multiplication with 10 ensures nonzero difference and should ensure stable computation
          value is not important for AMICI_ONE_STEP mode, only direction w.r.t. current t
          */
-        solver->AMISolve(std::max(*t,1.0) * 10, x, &dx, t,
-                                  AMICI_ONE_STEP);
+        if (it<1) {
+            int status = CVode(newton_sim.get(), std::max(*t,1.0) * 10, x->getNVector(), t, AMICI_ONE_STEP);
+            if(status != CV_SUCCESS)
+                throw CvodeException(status,"Error when calling CVode during Newton preequilibration simulation");
+        } else {
+            solver->AMISolve(std::max(*t,1.0) * 10, x, &dx, t, AMICI_ONE_STEP);
+        }
+
         model->fxdot(*t, x, &dx, &xdot);
         res_abs = sqrt(N_VDotProd(xdot.getNVector(), xdot.getNVector()));
         
@@ -289,9 +300,88 @@ void SteadystateProblem::getNewtonSimulation(ReturnData *rdata, Solver *solver,
         it_newton++;
         if (it_newton >= solver->getMaxSteps())
             throw NewtonFailure("Simulation based steady state failed to converge");
+    }
+}
+
+std::unique_ptr<void, std::function<void (void *)> > SteadystateProblem::createSteadystateSimSolver(
+        Solver *solver, Model *model, realtype tstart)
+{
+    /**
+     * New CVode object for preequilibration simulation is created
+     *
+     * @param[in] solver pointer to the AMICI solver object @type Solver
+     * @param[in] model pointer to the AMICI model object @type Model
+     * @param[in] tstart time point for starting Newton simulation @type realtype
+     */
+    
+    /* Create new CVode object */
+    auto newton_sim = std::unique_ptr<void, std::function<void(void *)>>(
+                CVodeCreate(solver->getLinearMultistepMethod(),
+                            solver->getNonlinearSolverIteration()),
+                [](void *ptr) { CVodeFree(&ptr); }
+    );
+    if (newton_sim == NULL)
+        throw AmiException("Failed to allocated solver memory!");
+    
+    int status = CVodeInit(newton_sim.get(), CVodeSolver::fxdot, RCONST(tstart), x->getNVector());
+    if(status != CV_SUCCESS)
+        throw CvodeException(status,"CVodeInit");
+    
+    /* Specify integration tolerances */
+    status = CVodeSStolerances(newton_sim.get(), RCONST(solver->getRelativeTolerance()),
+                               RCONST(solver->getAbsoluteTolerance()));
+    if(status != CV_SUCCESS)
+        throw CvodeException(status,"CVodeSStolerances");
+    
+    /* attaches userdata*/
+    status = CVodeSetUserData(newton_sim.get(), model);
+    if(status != CV_SUCCESS)
+        throw CvodeException(status,"CVodeSetUserData");
+    
+    /* specify maximal number of steps */
+    status = CVodeSetMaxNumSteps(newton_sim.get(), solver->getMaxSteps());
+    if(status != CV_SUCCESS)
+        throw CvodeException(status,"CVodeSetMaxNumSteps");
+    
+    /* activates stability limit detection */
+    status = CVodeSetStabLimDet(newton_sim.get(), solver->getStabilityLimitFlag());
+    if(status != CV_SUCCESS)
+        throw CvodeException(status,"CVodeSetStabLimDet");
+    
+    switch(solver->getLinearSolver()) {
             
+    case AMICI_DENSE:
+        /* Set up dense solver */
+        status = CVDense(newton_sim.get(), model->nx);
+        if(status != CV_SUCCESS)
+            throw CvodeException(status,"CVDense");
+
+        status = CVDlsSetDenseJacFn(newton_sim.get(), CVodeSolver::fJ);
+        if(status != CV_SUCCESS)
+            throw CvodeException(status,"CVDlsSetDenseJacFn");
+        break;
         
+    case AMICI_KLU:
+        /* Set up KLU solver */
+        status = CVKLU(newton_sim.get(), model->nx, model->nnz, CSC_MAT);
+        if(status != CV_SUCCESS)
+            throw CvodeException(status,"CVKLU");
+    
+        /* Pass Jacobian function to KLU solver */
+        status = CVSlsSetSparseJacFn(newton_sim.get(), CVodeSolver::fJSparse);
+        if(status != CV_SUCCESS)
+            throw CvodeException(status,"CVSlsSetSparseJacFn");
+    
+        /* Provide ordering to KLU solver */
+        status = CVKLUSetOrdering(newton_sim.get(), solver->getStateOrdering());
+        if(status != CV_SUCCESS)
+            throw CvodeException(status,"CVKLUSetOrdering");
+        break;
+            
+    default:
+        throw NewtonFailure("Invalid Choice of Solver!");
     }
+    return newton_sim;
 }
 
 } // namespace amici
diff --git a/tests/cpputest/expectedResults.h5 b/tests/cpputest/expectedResults.h5
index 1f5b32c772..07876b295f 100644
Binary files a/tests/cpputest/expectedResults.h5 and b/tests/cpputest/expectedResults.h5 differ
diff --git a/tests/cpputest/steadystate/tests1.cpp b/tests/cpputest/steadystate/tests1.cpp
index 2bdaf4b348..ca46fc4354 100644
--- a/tests/cpputest/steadystate/tests1.cpp
+++ b/tests/cpputest/steadystate/tests1.cpp
@@ -78,7 +78,6 @@ TEST(groupSteadystate, testReuseSolver) {
 }
 
 
-
 TEST(groupSteadystate, testSimulation) {
     amici::simulateVerifyWrite("/model_steadystate/nosensi/");
 }
@@ -109,4 +108,26 @@ TEST(groupSteadystate, testSensitivityForwardSPBCG) {
     amici::simulateVerifyWrite("/model_steadystate/nosensiSPBCG/", 10*TEST_ATOL, 10*TEST_RTOL);
 }
 
+TEST(groupSteadystate, testSensiFwdNewtonPreeq) {
+    amici::simulateVerifyWrite("/model_steadystate/sensifwdnewtonpreeq/");
+}
+
+TEST(groupSteadystate, testSensiAdjNewtonPreeq) {
+    amici::simulateVerifyWrite("/model_steadystate/sensiadjnewtonpreeq/");
+}
+
+TEST(groupSteadystate, testSensiFwdSimPreeq) {
+    amici::simulateVerifyWrite("/model_steadystate/sensifwdsimpreeq/");
+}
 
+TEST(groupSteadystate, testSensiAdjSimPreeq) {
+    amici::simulateVerifyWrite("/model_steadystate/sensiadjsimpreeq/");
+}
+
+TEST(groupSteadystate, testSensiFwdByhandPreeq) {
+    amici::simulateVerifyWrite("/model_steadystate/sensifwdbyhandpreeq/");
+}
+
+TEST(groupSteadystate, testSensiAdjByhandPreeq) {
+    amici::simulateVerifyWrite("/model_steadystate/sensiadjbyhandpreeq/");
+}
diff --git a/tests/cpputest/testOptions.h5 b/tests/cpputest/testOptions.h5
index 9d7d1ed038..6b782ea885 100644
Binary files a/tests/cpputest/testOptions.h5 and b/tests/cpputest/testOptions.h5 differ
diff --git a/tests/example_jakstat.py b/tests/example_jakstat.py
index 472bb686a0..2eca23cf1b 100755
--- a/tests/example_jakstat.py
+++ b/tests/example_jakstat.py
@@ -18,7 +18,7 @@ def __init__(self):
         self.numK = 2
 
         curPath = os.path.dirname(os.path.realpath(__file__))
-        dataPath = curPath + "/../examples/example_jakstat_adjoint/pnas_data_original.xls"
+        dataPath = curPath + "/../matlab/examples/example_jakstat_adjoint/pnas_data_original.xls"
         xls = pd.ExcelFile(dataPath).parse()
         self.modelOptions['ts'] = xls.time
         self.modelOptions['theta'] = np.array([0.60, 3, -0.95, -0.0075, 0,
diff --git a/tests/example_steadystate.py b/tests/example_steadystate.py
index 9374c684d5..b883e3766f 100755
--- a/tests/example_steadystate.py
+++ b/tests/example_steadystate.py
@@ -90,6 +90,210 @@ def writeSensiForwardErrorNewt(filename):
     ex.writeToFile(filename, '/model_steadystate/sensiforwarderrornewt/')
 
 
+def writeSensiFwdNewtonPreeq(filename):
+    ex = ExampleSteadystate()
+
+    ex.data = {}
+    ex.data['Y'] = np.array(
+        [[0.8410,1.6662,0.3813],
+        [0.7834,1.6230,0.2233],
+        [0.8187,1.5386,0.2208],
+        [0.7906,1.4476,0.2358],
+        [0.7184,1.3543,0.1409],
+        [0.6627,1.2840,0.1268],
+        [0.7099,1.1786,0.1289],
+        [0.7104,1.1362,0.1768],
+        [0.7089,1.0326,0.1127],
+        [0.6035,0.9752,0.0923]])
+    ex.data['Sigma_Y'] = np.ones(shape=ex.data['Y'].shape)
+    ex.data['Sigma_Z'] = []
+    ex.data['Z'] = []
+    ex.data['condition'] = ex.modelOptions['kappa']
+    ex.data['conditionPreequilibration'] = [0.1, 0.4, 0.7, 0.5]
+    ex.data['t'] = np.linspace(0, 5, 10)
+
+    ex.modelOptions['ts'] = np.linspace(0, 5, 10)
+    ex.solverOptions['sensi'] = 1
+    ex.solverOptions['sensi_meth'] = 1
+    ex.solverOptions['newton_preeq'] = True
+    ex.solverOptions['atol'] = 10**-16
+    ex.solverOptions['rtol'] = 10**-12
+    ex.solverOptions['sens_ind'] = np.arange(0, ex.numP)
+    ex.solverOptions['linsol'] = 9
+    ex.solverOptions['maxsteps'] = 10000
+    ex.solverOptions['newton_maxsteps'] = 20
+
+    ex.writeToFile(filename, '/model_steadystate/sensifwdnewtonpreeq/')
+
+def writeSensiAdjNewtonPreeq(filename):
+    ex = ExampleSteadystate()
+
+    ex.data = {}
+    ex.data['Y'] = np.array(
+        [[0.8410,1.6662,0.3813],
+        [0.7834,1.6230,0.2233],
+        [0.8187,1.5386,0.2208],
+        [0.7906,1.4476,0.2358],
+        [0.7184,1.3543,0.1409],
+        [0.6627,1.2840,0.1268],
+        [0.7099,1.1786,0.1289],
+        [0.7104,1.1362,0.1768],
+        [0.7089,1.0326,0.1127],
+        [0.6035,0.9752,0.0923]])
+    ex.data['Sigma_Y'] = np.ones(shape=ex.data['Y'].shape)
+    ex.data['Sigma_Z'] = []
+    ex.data['Z'] = []
+    ex.data['condition'] = ex.modelOptions['kappa']
+    ex.data['conditionPreequilibration'] = [0.1, 0.4, 0.7, 0.5]
+    ex.data['t'] = np.linspace(0, 5, 10)
+
+    ex.modelOptions['ts'] = np.linspace(0, 5, 10)
+    ex.solverOptions['sensi'] = 1
+    ex.solverOptions['sensi_meth'] = 2
+    ex.solverOptions['newton_preeq'] = True
+    ex.solverOptions['atol'] = 10**-16
+    ex.solverOptions['rtol'] = 10**-12
+    ex.solverOptions['sens_ind'] = np.arange(0, ex.numP)
+    ex.solverOptions['linsol'] = 9
+    ex.solverOptions['maxsteps'] = 10000
+    ex.solverOptions['newton_maxsteps'] = 20
+
+    ex.writeToFile(filename, '/model_steadystate/sensiadjnewtonpreeq/')
+
+
+def writeSensiFwdSimPreeq(filename):
+    ex = ExampleSteadystate()
+
+    ex.data = {}
+    ex.data['Y'] = np.array(
+        [[0.8410,1.6662,0.3813],
+        [0.7834,1.6230,0.2233],
+        [0.8187,1.5386,0.2208],
+        [0.7906,1.4476,0.2358],
+        [0.7184,1.3543,0.1409],
+        [0.6627,1.2840,0.1268],
+        [0.7099,1.1786,0.1289],
+        [0.7104,1.1362,0.1768],
+        [0.7089,1.0326,0.1127],
+        [0.6035,0.9752,0.0923]])
+    ex.data['Sigma_Y'] = np.ones(shape=ex.data['Y'].shape)
+    ex.data['Sigma_Z'] = []
+    ex.data['Z'] = []
+    ex.data['condition'] = ex.modelOptions['kappa']
+    ex.data['conditionPreequilibration'] = [0.1, 0.4, 0.7, 0.5]
+    ex.data['t'] = np.linspace(0, 5, 10)
+
+    ex.modelOptions['ts'] = np.linspace(0, 5, 10)
+    ex.solverOptions['sensi'] = 1
+    ex.solverOptions['sensi_meth'] = 1
+    ex.solverOptions['newton_preeq'] = True
+    ex.solverOptions['atol'] = 10 ** -16
+    ex.solverOptions['rtol'] = 10 ** -12
+    ex.solverOptions['sens_ind'] = np.arange(0, ex.numP)
+    ex.solverOptions['linsol'] = 9
+    ex.solverOptions['maxsteps'] = 10000
+    ex.solverOptions['newton_maxsteps'] = 0
+
+    ex.writeToFile(filename, '/model_steadystate/sensifwdsimpreeq/')
+
+
+def writeSensiAdjSimPreeq(filename):
+    ex = ExampleSteadystate()
+
+    ex.data = {}
+    ex.data['Y'] = np.array(
+        [[0.8410,1.6662,0.3813],
+        [0.7834,1.6230,0.2233],
+        [0.8187,1.5386,0.2208],
+        [0.7906,1.4476,0.2358],
+        [0.7184,1.3543,0.1409],
+        [0.6627,1.2840,0.1268],
+        [0.7099,1.1786,0.1289],
+        [0.7104,1.1362,0.1768],
+        [0.7089,1.0326,0.1127],
+        [0.6035,0.9752,0.0923]])
+    ex.data['Sigma_Y'] = np.ones(shape=ex.data['Y'].shape)
+    ex.data['Sigma_Z'] = []
+    ex.data['Z'] = []
+    ex.data['condition'] = ex.modelOptions['kappa']
+    ex.data['conditionPreequilibration'] = [0.1, 0.4, 0.7, 0.5]
+    ex.data['t'] = np.linspace(0, 5, 10)
+
+    ex.modelOptions['ts'] = np.linspace(0, 5, 10)
+    ex.solverOptions['sensi'] = 1
+    ex.solverOptions['sensi_meth'] = 2
+    ex.solverOptions['newton_preeq'] = True
+    ex.solverOptions['atol'] = 10 ** -16
+    ex.solverOptions['rtol'] = 10 ** -12
+    ex.solverOptions['sens_ind'] = np.arange(0, ex.numP)
+    ex.solverOptions['linsol'] = 9
+    ex.solverOptions['maxsteps'] = 10000
+    ex.solverOptions['newton_maxsteps'] = 0
+
+    ex.writeToFile(filename, '/model_steadystate/sensiadjsimpreeq/')
+
+def writeSensiFwdByhandPreeq(filename):
+    ex = ExampleSteadystate()
+
+    ex.modelOptions['ts'] = np.linspace(0, 5, 10)
+    ex.solverOptions['sensi'] = 1
+    ex.solverOptions['sensi_meth'] = 1
+    ex.solverOptions['newton_preeq'] = False
+    ex.solverOptions['x0'] = np.array([0.532609637980272, 0.625849232840357, 0.066666666666667])
+    ex.solverOptions['sx0'] = np.array(
+        [[-0.425457638009506, 0.499939012301881, 0], [-0.375463736779318, -0.999878024603762, -0.000000000000000],
+         [0.375463736779318, -0.441193089396203, 0], [0.300370989423454, 0.799902419683010, 0.000000000000000],
+         [0.425457638009506, 0.941132101698086, 0.153505672866270]])
+    ex.solverOptions['atol'] = 10 ** -16
+    ex.solverOptions['rtol'] = 10 ** -12
+    ex.solverOptions['sens_ind'] = np.arange(0, ex.numP)
+    ex.solverOptions['linsol'] = 9
+    ex.solverOptions['maxsteps'] = 10000
+    ex.solverOptions['newton_maxsteps'] = 0
+
+    ex.writeToFile(filename, '/model_steadystate/sensifwdbyhandpreeq/')
+
+
+def writeSensiAdjByhandPreeq(filename):
+    ex = ExampleSteadystate()
+
+    ex.data = {}
+    ex.data['Y'] = np.array(
+        [[0.8410,1.6662,0.3813],
+        [0.7834,1.6230,0.2233],
+        [0.8187,1.5386,0.2208],
+        [0.7906,1.4476,0.2358],
+        [0.7184,1.3543,0.1409],
+        [0.6627,1.2840,0.1268],
+        [0.7099,1.1786,0.1289],
+        [0.7104,1.1362,0.1768],
+        [0.7089,1.0326,0.1127],
+        [0.6035,0.9752,0.0923]])
+    ex.data['Sigma_Y'] = np.ones(shape=ex.data['Y'].shape)
+    ex.data['Sigma_Z'] = []
+    ex.data['Z'] = []
+    ex.data['condition'] = ex.modelOptions['kappa']
+    ex.data['t'] = np.linspace(0, 5, 10)
+
+    ex.modelOptions['ts'] = np.linspace(0, 5, 10)
+    ex.solverOptions['sensi'] = 1
+    ex.solverOptions['sensi_meth'] = 2
+    ex.solverOptions['newton_preeq'] = False
+    ex.solverOptions['x0'] = np.array([0.532609637980272, 0.625849232840357, 0.066666666666667])
+    ex.solverOptions['sx0'] = np.array(
+        [[-0.425457638009506, 0.499939012301881, 0], [-0.375463736779318, -0.999878024603762, -0.000000000000000],
+         [0.375463736779318, -0.441193089396203, 0], [0.300370989423454, 0.799902419683010, 0.000000000000000],
+         [0.425457638009506, 0.941132101698086, 0.153505672866270]])
+    ex.solverOptions['atol'] = 10 ** -16
+    ex.solverOptions['rtol'] = 10 ** -12
+    ex.solverOptions['sens_ind'] = np.arange(0, ex.numP)
+    ex.solverOptions['linsol'] = 9
+    ex.solverOptions['maxsteps'] = 10000
+    ex.solverOptions['newton_maxsteps'] = 0
+
+    ex.writeToFile(filename, '/model_steadystate/sensiadjbyhandpreeq/')
+
+
 def main():
     if len(sys.argv) < 2:
         print("Error: Must provide output file as first and only argument.")
@@ -103,6 +307,12 @@ def main():
     writeNosensiSPBCG(filename)
     writeSensiForwardErrorInt(filename)
     writeSensiForwardErrorNewt(filename)
+    writeSensiFwdNewtonPreeq(filename)
+    writeSensiAdjNewtonPreeq(filename)
+    writeSensiFwdSimPreeq(filename)
+    writeSensiAdjSimPreeq(filename)
+    writeSensiFwdByhandPreeq(filename)
+    writeSensiAdjByhandPreeq(filename)
     
 if __name__ == "__main__":
     main()