diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
index 01a0db41164cf071379df5caeb22197eedc9847a..411f9050013cf614d1b24c51f910be8eb24db30e 100644
--- a/src/finite-strain-elasticity.cc
+++ b/src/finite-strain-elasticity.cc
@@ -160,21 +160,16 @@ int main (int argc, char *argv[]) try
   // Make Python function that computes which vertices are on the Dirichlet boundary,
   // based on the vertex positions.
   std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+  auto pythonDirichletVertices = Python::make_function<bool>(Python::evaluate(lambda));
 
   // Same for the Neumann boundary
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
+  auto pythonNeumannVertices = Python::make_function<bool>(Python::evaluate(lambda));
 
   for (auto&& v : vertices(gridView))
   {
-    bool isDirichlet;
-    pythonDirichletVertices.evaluate(v.geometry().corner(0), isDirichlet);
-    dirichletVertices[indexSet.index(v)] = isDirichlet;
-
-    bool isNeumann;
-    pythonNeumannVertices.evaluate(v.geometry().corner(0), isNeumann);
-    neumannVertices[indexSet.index(v)] = isNeumann;
+    dirichletVertices[indexSet.index(v)] = pythonDirichletVertices(v.geometry().corner(0));
+    neumannVertices[indexSet.index(v)] = pythonNeumannVertices(v.geometry().corner(0));
   }
 
   BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
@@ -202,7 +197,7 @@ int main (int argc, char *argv[]) try
   SolutionType x(basis.size());
 
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > pythonInitialDeformation(Python::evaluate(lambda));
+  auto pythonInitialDeformation = Python::make_function<FieldVector<double,dim> >(Python::evaluate(lambda));
 
   Dune::Functions::interpolate(basis, x, pythonInitialDeformation);
 
@@ -328,7 +323,7 @@ int main (int argc, char *argv[]) try
     Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
 
     // Extract object member functions as Dune functions
-    PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> >   dirichletValues(dirichletValuesPythonObject.get("dirichletValues"));
+    auto dirichletValues = Python::make_function<FieldVector<double, dim> >(dirichletValuesPythonObject.get("dirichletValues"));
 
     Dune::Functions::interpolate(basis, x, dirichletValues, dirichletDofs);