diff --git a/include/amici/steadystateproblem.h b/include/amici/steadystateproblem.h
index f120094aa1..da6af249ac 100644
--- a/include/amici/steadystateproblem.h
+++ b/include/amici/steadystateproblem.h
@@ -378,15 +378,6 @@ class SteadystateProblem {
      */
     void getNewtonStep(Model& model);
 
-    /**
-     * @brief SUNDIALS context.
-     *
-     * We use a different context here, because we might
-     * use a different solver than the one used for the simulation, and we
-     * shouldn't use multiple solvers on the same context.
-     */
-    // TODO: attach error handler?
-    sundials::Context sunctx_;
     /** newton step */
     AmiVector delta_;
     /** previous newton step */
diff --git a/src/newton_solver.cpp b/src/newton_solver.cpp
index 5cad50d951..3238aea6fe 100644
--- a/src/newton_solver.cpp
+++ b/src/newton_solver.cpp
@@ -26,7 +26,9 @@ NewtonSolver::getSolver(Solver const& simulationSolver, Model const& model) {
 
     /* DIRECT SOLVERS */
     case LinearSolver::dense:
-        solver.reset(new NewtonSolverDense(model, simulationSolver.getSunContext()));
+        solver.reset(
+            new NewtonSolverDense(model, simulationSolver.getSunContext())
+        );
         break;
 
     case LinearSolver::band:
@@ -55,7 +57,9 @@ NewtonSolver::getSolver(Solver const& simulationSolver, Model const& model) {
     case LinearSolver::SuperLUMT:
         throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver");
     case LinearSolver::KLU:
-        solver.reset(new NewtonSolverSparse(model, simulationSolver.getSunContext()));
+        solver.reset(
+            new NewtonSolverSparse(model, simulationSolver.getSunContext())
+        );
         break;
     default:
         throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver");
diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp
index bf1c54ad4e..fe29ef4331 100644
--- a/src/steadystateproblem.cpp
+++ b/src/steadystateproblem.cpp
@@ -17,25 +17,29 @@ namespace amici {
 constexpr realtype conv_thresh = 1.0;
 
 SteadystateProblem::SteadystateProblem(Solver const& solver, Model const& model)
-    : delta_(model.nx_solver, sunctx_)
-    , delta_old_(model.nx_solver, sunctx_)
-    , ewt_(model.nx_solver, sunctx_)
-    , ewtQB_(model.nplist(), sunctx_)
-    , x_old_(model.nx_solver, sunctx_)
-    , xdot_(model.nx_solver, sunctx_)
-    , sdx_(model.nx_solver, model.nplist(), sunctx_)
-    , xB_(model.nJ * model.nx_solver, sunctx_)
-    , xQ_(model.nJ * model.nx_solver, sunctx_)
-    , xQB_(model.nplist(), sunctx_)
-    , xQBdot_(model.nplist(), sunctx_)
-    , steadystate_mask_(AmiVector(model.get_steadystate_mask(), sunctx_))
+    : delta_(model.nx_solver, solver.getSunContext())
+    , delta_old_(model.nx_solver, solver.getSunContext())
+    , ewt_(model.nx_solver, solver.getSunContext())
+    , ewtQB_(model.nplist(), solver.getSunContext())
+    , x_old_(model.nx_solver, solver.getSunContext())
+    , xdot_(model.nx_solver, solver.getSunContext())
+    , sdx_(model.nx_solver, model.nplist(), solver.getSunContext())
+    , xB_(model.nJ * model.nx_solver, solver.getSunContext())
+    , xQ_(model.nJ * model.nx_solver, solver.getSunContext())
+    , xQB_(model.nplist(), solver.getSunContext())
+    , xQBdot_(model.nplist(), solver.getSunContext())
+    , steadystate_mask_(
+          AmiVector(model.get_steadystate_mask(), solver.getSunContext())
+      )
     , max_steps_(solver.getNewtonMaxSteps())
     , dJydx_(model.nJ * model.nx_solver * model.nt(), 0.0)
     , state_(
-          {INFINITY,                                                 // t
-           AmiVector(model.nx_solver, sunctx_),                      // x
-           AmiVector(model.nx_solver, sunctx_),                      // dx
-           AmiVectorArray(model.nx_solver, model.nplist(), sunctx_), // sx
+          {INFINITY,                                           // t
+           AmiVector(model.nx_solver, solver.getSunContext()), // x
+           AmiVector(model.nx_solver, solver.getSunContext()), // dx
+           AmiVectorArray(
+               model.nx_solver, model.nplist(), solver.getSunContext()
+           ), // sx
            model.getModelState()}
       )
     , // state