diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index 6048c7d3d9ca019957375346b80e96df43fa3e93..6a18c1a51d16f9a4e379af97e5a03fb056295f63 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -1270,106 +1270,6 @@ assembleGradient(const std::vector<Configuration>& sol,
 
 }
 
-template <class GridType>
-void RodAssembler<GridType>::
-assembleGradientFD(const std::vector<Configuration>& sol,
-                 Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
-{
-    double eps = 1e-2;
-    if (grad.size() != sol.size())
-        grad.resize(sol.size());
-
-    // ///////////////////////////////////////////////////////////
-    //   Compute gradient by finite-difference approximation
-    // ///////////////////////////////////////////////////////////
-    std::vector<Configuration> forwardSolution = sol;
-    std::vector<Configuration> backwardSolution = sol;
-
-    // ////////////////////////////////////////////////////
-    //   Create local assembler
-    // ////////////////////////////////////////////////////
-
-    Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
-    Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
-    RodLocalStiffness<GridType,double> localStiffness(K, A);
-
-    // //////////////////////////////////////////////////////////////////
-    //   Store pointers to all elements so we can access them by index
-    // //////////////////////////////////////////////////////////////////
-
-    const typename GridType::Traits::LeafIndexSet& indexSet = grid_->leafIndexSet();
-
-    std::vector<EntityPointer> elements(grid_->size(0),grid_->template leafbegin<0>());
-    
-    typename GridType::template Codim<0>::LeafIterator eIt    = grid_->template leafbegin<0>();
-    typename GridType::template Codim<0>::LeafIterator eEndIt = grid_->template leafend<0>();
-
-    for (; eIt!=eEndIt; ++eIt)
-        elements[indexSet.index(*eIt)] = eIt;
-
-    Dune::array<Configuration,2> localReferenceConfiguration;
-    Dune::array<Configuration,2> localSolution;
-
-    // ///////////////////////////////////////////////////////////////
-    //   Loop over all blocks of the outer matrix
-    // ///////////////////////////////////////////////////////////////
-    for (size_t i=0; i<sol.size(); i++) {
-
-        for (int j=0; j<6; j++) {
-
-            double forwardEnergy = 0;
-            double backwardEnergy = 0;
-            
-            infinitesimalVariation(forwardSolution[i], eps, j);
-            infinitesimalVariation(backwardSolution[i], -eps, j);
-            
-            if (i != grid_->size(1)-1) {
-                
-                localReferenceConfiguration[0] = referenceConfiguration_[i];
-                localReferenceConfiguration[1] = referenceConfiguration_[i+1];
-                
-                localSolution[0] = forwardSolution[i];
-                localSolution[1] = forwardSolution[i+1];
-                
-                forwardEnergy += localStiffness.energy(*elements[i], localSolution, localReferenceConfiguration);
-                
-                localSolution[0] = backwardSolution[i];
-                localSolution[1] = backwardSolution[i+1];
-                
-                backwardEnergy += localStiffness.energy(*elements[i], localSolution, localReferenceConfiguration);
-                
-            } 
-            
-            if (i != 0) {
-                
-                localReferenceConfiguration[0] = referenceConfiguration_[i-1];
-                localReferenceConfiguration[1] = referenceConfiguration_[i];
-                
-                localSolution[0] = forwardSolution[i-1];
-                localSolution[1] = forwardSolution[i];
-                
-                forwardEnergy += localStiffness.energy(*elements[i-1], localSolution, localReferenceConfiguration);
-                
-                localSolution[0] = backwardSolution[i-1];
-                localSolution[1] = backwardSolution[i];
-                
-                backwardEnergy += localStiffness.energy(*elements[i-1], localSolution, localReferenceConfiguration);
-                
-            } 
-            
-            // Second derivative
-            grad[i][j] = (forwardEnergy - backwardEnergy) / (2*eps);
-            
-            forwardSolution[i]  = sol[i];
-            backwardSolution[i] = sol[i];
-            
-        }
-
-    }
-
-
-}
-
 
 template <class GridType>
 double RodAssembler<GridType>::
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index b8735f4529c169edfd6b7b27ab95a0c26d88ed36..a6f64872078c462ed5c6115c97af70f10f8e9ea2 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -262,9 +262,6 @@ public:
         void assembleGradient(const std::vector<Configuration>& sol,
                               Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
 
-        void assembleGradientFD(const std::vector<Configuration>& sol,
-                                Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
-
         /** \brief Compute the energy of a deformation state */
         double computeEnergy(const std::vector<Configuration>& sol) const;