diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index 41e4bf4aaf3d0e010ba967ec229db2e838af28b2..efaaa19a23205d255cd81d0a0203f90b4e0e82d5 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -516,4 +516,15 @@ EigenvalueSolver::_findLastRoots(Polynomial<2> P, const double epsilon) const
   return roots;
 }
 
+void
+EigenvalueSolver::computeExpForTinyMatrix(const TinyMatrix<3>& A, TinyMatrix<3>& expA)
+{
+  auto [eigenvalues, eigenvectors] = findEigen(A);
+  TinyMatrix<3> D                  = zero;
+  for (size_t i = 0; i < 3; ++i) {
+    D(i, i) = std::exp(eigenvalues[i]);
+  }
+  expA = eigenvectors * D * transpose(eigenvectors);
+}
+
 EigenvalueSolver::EigenvalueSolver() {}
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
index 6c01068e0c04d7fd89e97d2d45d8117367045934..728cacdfa59843cffc6d3d45e2ed0f716e03722a 100644
--- a/src/algebra/EigenvalueSolver.hpp
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -100,6 +100,8 @@ class EigenvalueSolver
 
   TinyVector<2> _findLastRoots(Polynomial<2> P, const double epsilon) const;
 
+  void computeExpForTinyMatrix(const TinyMatrix<3>& A, TinyMatrix<3>& expA);
+
   EigenvalueSolver();
   ~EigenvalueSolver() = default;
 };
diff --git a/src/scheme/HyperplasticSolver.cpp b/src/scheme/HyperplasticSolver.cpp
index e3f06ca9c279a8a4d91ed0ab6733dff29691ae41..9e8a272782d708c10a5a33d7398fba5ef5599f1e 100644
--- a/src/scheme/HyperplasticSolver.cpp
+++ b/src/scheme/HyperplasticSolver.cpp
@@ -203,9 +203,12 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
   double
   _compute_chi(const R3x3 S, const double yield) const
   {
-    const R3x3 Id                     = identity;
-    const R3x3 Sd                     = S - 1. / 3 * trace(S) * Id;
-    auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestB);
+    const R3x3 Id = identity;
+    const R3x3 Sd = S - 1. / 3 * trace(S) * Id;
+    //    auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestB);
+    R3x3 B;
+    EigenvalueSolver{}.computeExpForTinyMatrix(TestA2, B);
+    std::cout << B << "\n";
     exit(0);
     double chi         = 0;
     const double normS = frobeniusNorm(Sd);
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index 5efaa921b891b064679a3166d2d2221666661235..d2316d8845db03772ce0c0d7d56ec8b60bc3b96a 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -166,6 +166,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
       TinyMatrix<3> TestD{1e10, 0, 0, 0, 1, 0, 0, 0, 1};
       TinyMatrix<3> TestE{3e-10, 2e-10, 4e-10, 2e-10, 0, 2e-10, 4e-10, 2e-10, 3e-10};
 
+      TinyMatrix<3> expA2;
+
       auto [eigenvalues, eigenmatrix] = EigenvalueSolver{}.findEigen(TestA);
       TinyMatrix<3> Diag              = zero;
       for (size_t i = 0; i < 3; ++i) {
@@ -177,7 +179,6 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
           REQUIRE((B(i, j) - TestA(i, j)) / frobeniusNorm(TestA) == Catch::Approx(0).margin(1E-8));
         }
       }
-
       auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestA2);
       Diag                              = zero;
       for (size_t i = 0; i < 3; ++i) {
@@ -189,6 +190,17 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
           REQUIRE((B(i, j) - TestA2(i, j)) / frobeniusNorm(TestA2) == Catch::Approx(0).margin(1E-8));
         }
       }
+      EigenvalueSolver{}.computeExpForTinyMatrix(TestA2, expA2);
+      for (size_t i = 0; i < 3; ++i) {
+        Diag(i, i) = std::exp(eigenvalues2[i]);
+      }
+      B = eigenmatrix2 * Diag * transpose(eigenmatrix2);
+
+      for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; j < 3; ++j) {
+          REQUIRE(expA2(i, j) - B(i, j) == Catch::Approx(0).margin(1E-12));
+        }
+      }
 
       auto [eigenvalues3, eigenmatrix3] = EigenvalueSolver{}.findEigen(TestB);
       Diag                              = zero;