diff --git a/dune/elasticity/Makefile.am b/dune/elasticity/Makefile.am
index 19010a7..70e2d68 100644
--- a/dune/elasticity/Makefile.am
+++ b/dune/elasticity/Makefile.am
@@ -2,6 +2,7 @@ lib_LTLIBRARIES = libelasticityupscale.la
 
 libelasticityupscale_la_SOURCES = \
         boundarygrid.cpp \
+	dynmatrixev.cpp \
         material.cpp \
         materials.cpp \
         mpc.cpp
diff --git a/dune/elasticity/dynmatrixev.cpp b/dune/elasticity/dynmatrixev.cpp
new file mode 100644
index 0000000..d5880a9
--- /dev/null
+++ b/dune/elasticity/dynmatrixev.cpp
@@ -0,0 +1,152 @@
+#ifndef DUNE_FMATRIXEIGENVALUES_EXT_CC 
+#define DUNE_FMATRIXEIGENVALUES_EXT_CC 
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <cmath>
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+
+#if HAVE_LAPACK
+
+// nonsymetric matrices
+#define DGEEV_FORTRAN FC_FUNC (dgeev, DGEEV)
+
+// dsyev declaration (in liblapack)
+extern "C" {
+
+/*
+    *
+    **  purpose
+    **  =======
+    **
+    **  xgeev computes for an N-by-N BASE DATA TYPE nonsymmetric matrix A, the
+    **  eigenvalues and, optionally, the left and/or right eigenvectors.
+    **
+    **  The right eigenvector v(j) of A satisfies
+    **                   A * v(j) = lambda(j) * v(j)
+    **  where lambda(j) is its eigenvalue.
+    **  The left eigenvector u(j) of A satisfies
+    **                u(j)**T * A = lambda(j) * u(j)**T
+    **  where u(j)**T denotes the transpose of u(j).
+    **
+    **  The computed eigenvectors are normalized to have Euclidean norm
+    **  equal to 1 and largest component real.
+    **
+    **  arguments
+    **  =========
+    **
+    **  jobvl   (input) char
+    **          = 'n': left eigenvectors of a are not computed;
+    **          = 'v': left eigenvectors of a are computed.
+    **
+    **  jobvr   (input) char
+    **          = 'n': right eigenvectors of a are not computed;
+    **          = 'v': right eigenvectors of a are computed.
+    **
+    **  n       (input) long int
+    **          the order of the matrix v. v >= 0.
+    **
+    **  a       (input/output) BASE DATA TYPE array, dimension (lda,n)
+    **          on entry, the n-by-n matrix a.
+    **          on exit, a has been overwritten.
+    **
+    **  lda     (input) long int
+    **          the leading dimension of the array a.  lda >= max(1,n).
+    **
+    **  wr      (output) BASE DATA TYPE array, dimension (n)
+    **  wi      (output) BASE DATA TYPE array, dimension (n)
+    **          wr and wi contain the real and imaginary parts,
+    **          respectively, of the computed eigenvalues.  complex
+    **          conjugate pairs of eigenvalues appear consecutively
+    **          with the eigenvalue having the positive imaginary part
+    **          first.
+    **
+    **  vl      (output) COMPLEX DATA TYPE array, dimension (ldvl,n)
+    **          if jobvl = 'v', the left eigenvectors u(j) are stored one
+    **          after another in the columns of vl, in the same order
+    **          as their eigenvalues.
+    **          if jobvl = 'n', vl is not referenced.
+    **          if the j-th eigenvalue is real, then u(j) = vl(:,j),
+    **          the j-th column of vl.
+    **          if the j-th and (j+1)-st eigenvalues form a complex
+    **          conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and
+    **          u(j+1) = vl(:,j) - i*vl(:,j+1).
+    **
+    **  ldvl    (input) long int
+    **          the leading dimension of the array vl.  ldvl >= 1; if
+    **          jobvl = 'v', ldvl >= n.
+    **
+    **  vr      (output) COMPLEX DATA TYPE array, dimension (ldvr,n)
+    **          if jobvr = 'v', the right eigenvectors v(j) are stored one
+    **          after another in the columns of vr, in the same order
+    **          as their eigenvalues.
+    **          if jobvr = 'n', vr is not referenced.
+    **          if the j-th eigenvalue is real, then v(j) = vr(:,j),
+    **          the j-th column of vr.
+    **          if the j-th and (j+1)-st eigenvalues form a complex
+    **          conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and
+    **          v(j+1) = vr(:,j) - i*vr(:,j+1).
+    **
+    **  ldvr    (input) long int 
+    **          the leading dimension of the array vr.  ldvr >= 1; if
+    **          jobvr = 'v', ldvr >= n.
+    **
+    **  work    (workspace/output) BASE DATA TYPE array, dimension (max(1,lwork))
+    **          on exit, if info = 0, work(1) returns the optimal lwork.
+    **
+    **  lwork   (input) long int
+    **          the dimension of the array work.  lwork >= max(1,3*n), and
+    **          if jobvl = 'v' or jobvr = 'v', lwork >= 4*n.  for good
+    **          performance, lwork must generally be larger.
+    **
+    **          if lwork = -1, then a workspace query is assumed; the routine
+    **          only calculates the optimal size of the work array, returns
+    **          this value as the first entry of the work array, and no error
+    **          message related to lwork is issued by xerbla.
+    **
+    **  info    (output) long int
+    **          = 0:  successful exit
+    **          < 0:  if info = -i, the i-th argument had an illegal value.
+    **          > 0:  if info = i, the qr algorithm failed to compute all the
+    **                eigenvalues, and no eigenvectors have been computed;
+    **                elements i+1:n of wr and wi contain eigenvalues which
+    **                have converged.
+    **
+**/
+
+extern void DGEEV_FORTRAN(const char* jobvl, const char* jobvr, const long
+    int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+    const long int* ldvl, double* vr, const long int* ldvr, double* work,
+    const long int* lwork, const long int* info);
+
+} // end extern C 
+#endif
+
+namespace Dune {
+
+namespace DynamicMatrixHelp {
+
+void eigenValuesNonsymLapackCall(
+      const char* jobvl, const char* jobvr, const long
+      int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+      const long int* ldvl, double* vr, const long int* ldvr, double* work,
+      const long int* lwork, const long int* info)
+{
+#if HAVE_LAPACK 
+  // call LAPACK dgeev 
+  DGEEV_FORTRAN(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr,
+                work, lwork, info);
+#else 
+  DUNE_THROW(NotImplemented,"eigenValuesNonsymLapackCall: LAPACK not found!");
+#endif
+}
+
+} // end namespace FMatrixHelp 
+
+} // end namespace Dune 
+#endif
diff --git a/dune/elasticity/dynmatrixev.hh b/dune/elasticity/dynmatrixev.hh
new file mode 100644
index 0000000..570c5cc
--- /dev/null
+++ b/dune/elasticity/dynmatrixev.hh
@@ -0,0 +1,67 @@
+#ifndef DUNE_DYNMATRIXEIGENVALUES_HH
+#define DUNE_DYNMATRIXEIGENVALUES_HH
+
+#include <dune/common/dynmatrix.hh>
+
+namespace Dune {
+
+namespace DynamicMatrixHelp {
+
+// defined in fmatrixev_ext.cpp
+extern void eigenValuesNonsymLapackCall(
+    const char* jobvl, const char* jobvr, const long
+    int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+    const long int* ldvl, double* vr, const long int* ldvr, double* work,
+    const long int* lwork, const long int* info);
+
+template <typename K, class C>
+static void eigenValuesNonSym(const DynamicMatrix<K>& matrix,
+                              DynamicVector<C>& eigenValues)
+{
+  {
+    const long int N = matrix.rows();
+    const char jobvl = 'n';
+    const char jobvr = 'n';
+
+    // matrix to put into dgeev 
+    double matrixVector[N * N]; 
+
+    // copy matrix  
+    int row = 0;
+    for(int i=0; i<N; ++i) 
+    {
+      for(int j=0; j<N; ++j, ++row) 
+      {
+        matrixVector[ row ] = matrix[ i ][ j ];
+      }
+    }
+
+    // working memory 
+    double eigenR[N]; 
+    double eigenI[N];
+    double work[3*N];
+
+    // return value information 
+    long int info = 0;
+    long int lwork = 3*N;
+
+    // call LAPACK routine (see fmatrixev_ext.cc)
+    eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N, 
+                                &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
+                                &lwork, &info);
+
+    if( info != 0 ) 
+    {
+      std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
+      DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
+    }
+    for (int i=0;i<N;++i)
+      eigenValues[i] = std::complex<double>(eigenR[i], eigenI[i]);
+  }
+}
+
+}
+
+}
+
+#endif
diff --git a/dune/elasticity/fmatrixev_ext.cc b/dune/elasticity/fmatrixev_ext.cc
new file mode 100644
index 0000000..a80c060
--- /dev/null
+++ b/dune/elasticity/fmatrixev_ext.cc
@@ -0,0 +1,152 @@
+#ifndef DUNE_FMATRIXEIGENVALUES_EXT_CC 
+#define DUNE_FMATRIXEIGENVALUES_EXT_CC 
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <cmath>
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+
+#if HAVE_LAPACK
+
+// nonsymetric matrices
+#define DGEEV_FORTRAN FC_FUNC (dgeev, DGEEV)
+
+// dsyev declaration (in liblapack)
+extern "C" {
+
+/*
+    *
+    **  purpose
+    **  =======
+    **
+    **  xgeev computes for an N-by-N BASE DATA TYPE nonsymmetric matrix A, the
+    **  eigenvalues and, optionally, the left and/or right eigenvectors.
+    **
+    **  The right eigenvector v(j) of A satisfies
+    **                   A * v(j) = lambda(j) * v(j)
+    **  where lambda(j) is its eigenvalue.
+    **  The left eigenvector u(j) of A satisfies
+    **                u(j)**T * A = lambda(j) * u(j)**T
+    **  where u(j)**T denotes the transpose of u(j).
+    **
+    **  The computed eigenvectors are normalized to have Euclidean norm
+    **  equal to 1 and largest component real.
+    **
+    **  arguments
+    **  =========
+    **
+    **  jobvl   (input) char
+    **          = 'n': left eigenvectors of a are not computed;
+    **          = 'v': left eigenvectors of a are computed.
+    **
+    **  jobvr   (input) char
+    **          = 'n': right eigenvectors of a are not computed;
+    **          = 'v': right eigenvectors of a are computed.
+    **
+    **  n       (input) long int
+    **          the order of the matrix v. v >= 0.
+    **
+    **  a       (input/output) BASE DATA TYPE array, dimension (lda,n)
+    **          on entry, the n-by-n matrix a.
+    **          on exit, a has been overwritten.
+    **
+    **  lda     (input) long int
+    **          the leading dimension of the array a.  lda >= max(1,n).
+    **
+    **  wr      (output) BASE DATA TYPE array, dimension (n)
+    **  wi      (output) BASE DATA TYPE array, dimension (n)
+    **          wr and wi contain the real and imaginary parts,
+    **          respectively, of the computed eigenvalues.  complex
+    **          conjugate pairs of eigenvalues appear consecutively
+    **          with the eigenvalue having the positive imaginary part
+    **          first.
+    **
+    **  vl      (output) COMPLEX DATA TYPE array, dimension (ldvl,n)
+    **          if jobvl = 'v', the left eigenvectors u(j) are stored one
+    **          after another in the columns of vl, in the same order
+    **          as their eigenvalues.
+    **          if jobvl = 'n', vl is not referenced.
+    **          if the j-th eigenvalue is real, then u(j) = vl(:,j),
+    **          the j-th column of vl.
+    **          if the j-th and (j+1)-st eigenvalues form a complex
+    **          conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and
+    **          u(j+1) = vl(:,j) - i*vl(:,j+1).
+    **
+    **  ldvl    (input) long int
+    **          the leading dimension of the array vl.  ldvl >= 1; if
+    **          jobvl = 'v', ldvl >= n.
+    **
+    **  vr      (output) COMPLEX DATA TYPE array, dimension (ldvr,n)
+    **          if jobvr = 'v', the right eigenvectors v(j) are stored one
+    **          after another in the columns of vr, in the same order
+    **          as their eigenvalues.
+    **          if jobvr = 'n', vr is not referenced.
+    **          if the j-th eigenvalue is real, then v(j) = vr(:,j),
+    **          the j-th column of vr.
+    **          if the j-th and (j+1)-st eigenvalues form a complex
+    **          conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and
+    **          v(j+1) = vr(:,j) - i*vr(:,j+1).
+    **
+    **  ldvr    (input) long int 
+    **          the leading dimension of the array vr.  ldvr >= 1; if
+    **          jobvr = 'v', ldvr >= n.
+    **
+    **  work    (workspace/output) BASE DATA TYPE array, dimension (max(1,lwork))
+    **          on exit, if info = 0, work(1) returns the optimal lwork.
+    **
+    **  lwork   (input) long int
+    **          the dimension of the array work.  lwork >= max(1,3*n), and
+    **          if jobvl = 'v' or jobvr = 'v', lwork >= 4*n.  for good
+    **          performance, lwork must generally be larger.
+    **
+    **          if lwork = -1, then a workspace query is assumed; the routine
+    **          only calculates the optimal size of the work array, returns
+    **          this value as the first entry of the work array, and no error
+    **          message related to lwork is issued by xerbla.
+    **
+    **  info    (output) long int
+    **          = 0:  successful exit
+    **          < 0:  if info = -i, the i-th argument had an illegal value.
+    **          > 0:  if info = i, the qr algorithm failed to compute all the
+    **                eigenvalues, and no eigenvectors have been computed;
+    **                elements i+1:n of wr and wi contain eigenvalues which
+    **                have converged.
+    **
+**/
+
+extern void DGEEV_FORTRAN(const char* jobvl, const char* jobvr, const long
+    int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+    const long int* ldvl, double* vr, const long int* ldvr, double* work,
+    const long int* lwork, const long int* info);
+
+} // end extern C 
+#endif
+
+namespace Dune {
+
+namespace FMatrixHelp {
+
+void eigenValuesNonsymLapackCall(
+      const char* jobvl, const char* jobvr, const long
+      int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+      const long int* ldvl, double* vr, const long int* ldvr, double* work,
+      const long int* lwork, const long int* info)
+{
+#if HAVE_LAPACK 
+  // call LAPACK dgeev 
+  DGEEV_FORTRAN(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr,
+                work, lwork, info);
+#else 
+  DUNE_THROW(NotImplemented,"eigenValuesNonsymLapackCall: LAPACK not found!");
+#endif
+}
+
+} // end namespace FMatrixHelp 
+
+} // end namespace Dune 
+#endif
diff --git a/dune/elasticity/fmatrixev_ext.hh b/dune/elasticity/fmatrixev_ext.hh
new file mode 100644
index 0000000..bfdebcc
--- /dev/null
+++ b/dune/elasticity/fmatrixev_ext.hh
@@ -0,0 +1,70 @@
+#ifndef DUNE_FMATRIXEIGENVALUES_EXT_HH
+#define DUNE_FMATRIXEIGENVALUES_EXT_HH
+
+namespace Dune {
+
+namespace FMatrixHelp {
+
+// defined in fmatrixev_ext.cpp
+extern void eigenValuesNonsymLapackCall(
+    const char* jobvl, const char* jobvr, const long
+    int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+    const long int* ldvl, double* vr, const long int* ldvr, double* work,
+    const long int* lwork, const long int* info);
+
+template <int dim, typename K, class C>
+static void eigenValuesNonSym(const FieldMatrix<K, dim, dim>& matrix,
+                              FieldVector<C, dim>& eigenValues)
+{
+  {
+    const long int N = dim ;
+    const char jobvl = 'n';
+    const char jobvr = 'n';
+
+    // length of matrix vector 
+    const long int w = N * N ;
+
+    // matrix to put into dgeev 
+    double matrixVector[dim * dim]; 
+
+    // copy matrix  
+    int row = 0;
+    for(int i=0; i<dim; ++i) 
+    {
+      for(int j=0; j<dim; ++j, ++row) 
+      {
+        matrixVector[ row ] = matrix[ i ][ j ];
+      }
+    }
+
+    // working memory 
+    double eigenR[dim]; 
+    double eigenI[dim];
+    double work[3*dim];
+
+    // return value information 
+    long int info = 0;
+    long int lwork = 3*dim;
+
+    // call LAPACK routine (see fmatrixev_ext.cc)
+    eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N, 
+                                &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
+                                &lwork, &info);
+
+    if( info != 0 ) 
+    {
+      std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
+      DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
+    }
+    for (int i=0;i<N;++i) {
+      eigenValues[i].real = eigenR[i];
+      eigenValues[i].imag = eigenI[i];
+    }
+  }
+}
+
+}
+
+}
+
+#endif