diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 3d3fa4d8def6a22dd602b2bcc046a767a8f30463..6348e753bf53c8ba3469fcbfd18ef9461019d43f 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -60,6 +60,7 @@
 #include <dune/gfe/localgeodesicfeadolcstiffness.hh>
 #include <dune/gfe/cosseratvtkwriter.hh>
 #include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/riemannianpnsolver.hh>
 #include <dune/gfe/riemanniantrsolver.hh>
 #include <dune/gfe/vertexnormals.hh>
 #include <dune/gfe/surfacecosseratenergy.hh>
@@ -152,8 +153,9 @@ int main (int argc, char *argv[]) try
   int numLevels                         = parameterSet.get<int>("numLevels");
   int numHomotopySteps                  = parameterSet.get<int>("numHomotopySteps");
   const double tolerance                = parameterSet.get<double>("tolerance");
-  const int maxTrustRegionSteps         = parameterSet.get<int>("maxTrustRegionSteps");
+  const int maxSolverSteps              = parameterSet.get<int>("maxSolverSteps");
   const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+  const double initialRegularization    = parameterSet.get<double>("initialRegularization");
   const int multigridIterations         = parameterSet.get<int>("numIt");
   const int nu1                         = parameterSet.get<int>("nu1");
   const int nu2                         = parameterSet.get<int>("nu2");
@@ -577,26 +579,6 @@ int main (int argc, char *argv[]) try
 
     GeodesicFEAssembler<FEBasis,TargetSpace> assembler(gridView, &localGFEADOLCStiffness);
 
-    // /////////////////////////////////////////////////
-    //   Create a Riemannian trust-region solver
-    // /////////////////////////////////////////////////
-
-    RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
-    solver.setup(*grid,
-                 &assembler,
-                 x,
-                 dirichletDofs,
-                 tolerance,
-                 maxTrustRegionSteps,
-                 initialTrustRegionRadius,
-                 multigridIterations,
-                 mgTolerance,
-                 mu, nu1, nu2,
-                 baseIterations,
-                 baseTolerance,
-                 instrumented);
-
-    solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
 
     ////////////////////////////////////////////////////////
     //   Set Dirichlet values
@@ -625,16 +607,50 @@ int main (int argc, char *argv[]) try
         x[j].r = ddV[j];
         x[j].q.set(dOV[j]);
       }
+    // /////////////////////////////////////////////////
+    //   Create the solver and solve
+    // /////////////////////////////////////////////////
 
+    if (parameterSet.get<std::string>("solvertype") == "multigrid") {
+      RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
+      solver.setup(*grid,
+                   &assembler,
+                   x,
+                   dirichletDofs,
+                   tolerance,
+                   maxSolverSteps,
+                   initialTrustRegionRadius,
+                   multigridIterations,
+                   mgTolerance,
+                   mu, nu1, nu2,
+                   baseIterations,
+                   baseTolerance,
+                   instrumented);
+
+      solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+      solver.setInitialIterate(x);
+      solver.solve();
+
+      x = solver.getSol();
+    } else { //parameterSet.get<std::string>("solvertype") == "cholmod"
+      RiemannianProximalNewtonSolver<FEBasis,TargetSpace> solver;
+      solver.setup(*grid,
+                   &assembler,
+                   x,
+                   dirichletDofs,
+                   tolerance,
+                   maxSolverSteps,
+                   initialRegularization,
+                   instrumented);
+      solver.setInitialIterate(x);
+      solver.solve();
+
+      x = solver.getSol();
+    }
     // /////////////////////////////////////////////////////
     //   Solve!
     // /////////////////////////////////////////////////////
 
-    solver.setInitialIterate(x);
-    solver.solve();
-
-    x = solver.getSol();
-
     std::cout << "Overall calculation took " << overallTimer.elapsed() << " sec." << std::endl;
 
     /////////////////////////////////
diff --git a/src/film-on-substrate.parset b/src/film-on-substrate.parset
index 645c408bf9efca144da5cdc2c0dfe1b101901442..211041833d47d7f066c91dcc8105dbfc5ae427b6 100644
--- a/src/film-on-substrate.parset
+++ b/src/film-on-substrate.parset
@@ -50,14 +50,31 @@ initialDeformation = "x"
 #  Solver parameters
 #############################################
 
+# Inner solver, cholmod or multigrid
+solvertype = cholmod
+
 # Number of homotopy steps for the Dirichlet boundary conditions
 numHomotopySteps = 1
 
-# Tolerance of the trust region solver
+# Tolerance of the solver
 tolerance = 1e-3
 
-# Max number of steps of the trust region solver
-maxTrustRegionSteps = 20
+# Max number of solver steps
+maxSolverSteps = 300
+
+# Measure convergence
+instrumented = 0
+
+#############################################
+#  Solver parameters specific for proximal newton solver using cholmod
+#############################################
+
+# initial regularization parameter
+initialRegularization = 10000000000
+
+#############################################
+#  Solver parameters specific for trust-region solver using multigrid solver 
+#############################################
 
 trustRegionScaling = 1 1 1 0.01 0.01 0.01