diff --git a/src/timestepping/eulerpair.cc b/src/timestepping/eulerpair.cc
index 364e64fd3f3f2bfba510f8d8640bcb7e2db8d48b..dfef35b84305d13bff5f0897b6f19cbbb18640fd 100644
--- a/src/timestepping/eulerpair.cc
+++ b/src/timestepping/eulerpair.cc
@@ -1,12 +1,11 @@
-
 template <class VectorType, class MatrixType, class FunctionType, int dim>
 EulerPair<VectorType, MatrixType, FunctionType, dim>::EulerPair(
-    MatrixType const &_A, MatrixType const &_B, VectorType const &_u_initial,
+    MatrixType const &_A, MatrixType const &_M, VectorType const &_u_initial,
     VectorType const &_v_initial,
     Dune::BitSetVector<dim> const &_dirichletNodes,
     FunctionType const &_dirichletFunction)
     : A(_A),
-      B(_B),
+      M(_M),
       u(_u_initial),
       v(_v_initial),
       dirichletNodes(_dirichletNodes),
@@ -21,26 +20,56 @@ void EulerPair<VectorType, MatrixType, FunctionType, dim>::nextTimeStep() {
 template <class VectorType, class MatrixType, class FunctionType, int dim>
 void EulerPair<VectorType, MatrixType, FunctionType, dim>::setup(
     VectorType const &ell, double _tau, double time, VectorType &problem_rhs,
-    VectorType &problem_iterate, MatrixType &problem_AB) {
+    VectorType &problem_iterate, MatrixType &problem_AM) {
   postProcessCalled = false;
 
   tau = _tau;
 
-  problem_rhs = ell;
-  Arithmetic::addProduct(problem_rhs, 1.0 / tau, B, v_o);
-  Arithmetic::subtractProduct(problem_rhs, A, u_o);
+  MatrixType const &C = A;
+  double const wc = 0.0;
+  /* We start out with the formulation
+
+       M a + C v + A u = ell
+
+     Backward Euler means
+
+       a1 = 1.0/tau ( v1 - v0 )
+       u1 = tau v1 + u0
+
+     in summary, we get at time t=1
+
+       M [1.0/tau ( v1 - v0 )] + C v1
+       + A [tau v1 + u0] = ell
+
+     or
+
+       1.0/tau M v1 + C v1 + tau A v1
+       = [1.0/tau M + C + tau A] v1
+       = ell + 1.0/tau M v0 - A u0
+  */
+
+  // set up LHS (for fixed tau, we'd only really have to do this once)
+  {
+    Dune::MatrixIndexSet indices(A.N(), A.M());
+    indices.import(A);
+    indices.import(M);
+    indices.import(C);
+    indices.exportIdx(problem_AM);
+  }
+  problem_AM = 0.0;
+  Arithmetic::addProduct(problem_AM, (1.0 - wc) / tau, M);
+  Arithmetic::addProduct(problem_AM, wc, C);
+  Arithmetic::addProduct(problem_AM, tau, A);
 
-  // For fixed tau, we'd only really have to do this once
-  Dune::MatrixIndexSet indices(A.N(), A.M());
-  indices.import(A);
-  indices.import(B);
-  indices.exportIdx(problem_AB);
-  problem_AB = 0.0;
-  Arithmetic::addProduct(problem_AB, tau, A);
-  Arithmetic::addProduct(problem_AB, 1.0 / tau, B);
+  // set up RHS
+  {
+    problem_rhs = ell;
+    Arithmetic::addProduct(problem_rhs, (1.0 - wc) / tau, M, v_o);
+    Arithmetic::subtractProduct(problem_rhs, A, u_o);
+  }
 
   // v_o makes a good initial iterate; we could use anything, though
-  problem_iterate = v_o;
+  problem_iterate = 0.0;
 
   for (size_t i = 0; i < dirichletNodes.size(); ++i)
     switch (dirichletNodes[i].count()) {
diff --git a/src/timestepping/eulerpair.hh b/src/timestepping/eulerpair.hh
index 289f98d26f1278d36e91d32047257e465999cb3e..105529e19971fef6f5e3284570a740d1e8aebd2a 100644
--- a/src/timestepping/eulerpair.hh
+++ b/src/timestepping/eulerpair.hh
@@ -5,7 +5,7 @@ template <class VectorType, class MatrixType, class FunctionType, int dim>
 class EulerPair
     : public TimeSteppingScheme<VectorType, MatrixType, FunctionType, dim> {
 public:
-  EulerPair(MatrixType const &_A, MatrixType const &_B,
+  EulerPair(MatrixType const &_A, MatrixType const &_M,
             VectorType const &_u_initial, VectorType const &_v_initial,
             Dune::BitSetVector<dim> const &_dirichletNodes,
             FunctionType const &_dirichletFunction);
@@ -19,7 +19,7 @@ class EulerPair
 
 private:
   MatrixType const &A;
-  MatrixType const &B;
+  MatrixType const &M;
   VectorType u;
   VectorType v;
   Dune::BitSetVector<dim> const &dirichletNodes;