diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index b232e988f186054f5ba9990ac23a36d737267261..02d2d2793805fe0530ee0a113ff6fed10d05723f 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -78,6 +78,13 @@ public:
         GlobalMatrix& Jac, GlobalVector& res,
         GlobalVector& minus_delta_x) const = 0;
 
+    //! Apply known solutions to the linearized equation system
+    //! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
+    //! \pre computeKnownSolutions() must have been called before.
+    virtual void applyKnownSolutionsPETScSNES(GlobalMatrix& Jac,
+                                              GlobalVector& res,
+                                              GlobalVector& x) const = 0;
+
     virtual void updateConstraints(GlobalVector& /*lower*/,
                                    GlobalVector& /*upper*/,
                                    int /*process_id*/) = 0;
diff --git a/NumLib/ODESolver/PETScNonlinearSolver.cpp b/NumLib/ODESolver/PETScNonlinearSolver.cpp
index 7eb3475b026ad2da2338af942f9f7daddb8ebf29..193b0dc65367b0a60f66f189ddfc4af8ac802e07 100644
--- a/NumLib/ODESolver/PETScNonlinearSolver.cpp
+++ b/NumLib/ODESolver/PETScNonlinearSolver.cpp
@@ -67,7 +67,7 @@ PetscErrorCode updateResidual(SNES /*snes*/, Vec x, Vec petsc_r,
     context->J->finalizeAssembly();
 
     context->system->getJacobian(*context->J);
-    context->system->applyKnownSolutionsNewton(
+    context->system->applyKnownSolutionsPETScSNES(
         *context->J, *context->r, *context->x[context->process_id]);
 
     VecCopy(context->r->getRawVector(), petsc_r);
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index d7e531afdc305b239d9b8a8f7696eca768cb0e94..e618773be5c518e8467a63d8556474d79aabe8fa 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -152,6 +152,28 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     MathLib::applyKnownSolution(Jac, res, minus_delta_x, ids, values);
 }
 
+void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
+                              NonlinearSolverTag::Newton>::
+    applyKnownSolutionsPETScSNES(GlobalMatrix& Jac, GlobalVector& res,
+                                 GlobalVector& x) const
+{
+    if (!_known_solutions)
+    {
+        return;
+    }
+
+    using IndexType = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
+    std::vector<IndexType> ids;
+    for (auto const& bc : *_known_solutions)
+    {
+        ids.insert(end(ids), begin(bc.ids), end(bc.ids));
+    }
+
+    // For the Newton method the values must be zero
+    std::vector<double> values(ids.size(), 0);
+    MathLib::applyKnownSolution(Jac, res, x, ids, values);
+}
+
 TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
                          NonlinearSolverTag::Picard>::
     TimeDiscretizedODESystem(const int process_id, ODE& ode,
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index b53f6ca426bebc5c8e6390264fd460e14573a6c8..46dd959689fa3d532fcb08c3904ec13db048179b 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -108,6 +108,9 @@ public:
     void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
                                    GlobalVector& minus_delta_x) const override;
 
+    void applyKnownSolutionsPETScSNES(GlobalMatrix& Jac, GlobalVector& res,
+                                      GlobalVector& x) const override;
+
     void updateConstraints(GlobalVector& lower, GlobalVector& upper,
                            int const process_id) override
     {