diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index c71be58ddcb1c37498a172f30a1eb39b89f43aab..842150bb5e2eabb17da38c0238cb57499f631c51 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -1,9 +1,9 @@
 #include <algebra/EigenvalueSolver.hpp>
+#include <cmath>
 #include <utils/pugs_config.hpp>
-
 #ifdef PUGS_HAS_SLEPC
 #include <slepc.h>
-
+#include <slepceps.h>
 struct EigenvalueSolver::Internals
 {
   static PetscReal
@@ -152,6 +152,43 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
   EPSDestroy(&eps);
 }
 
+void
+EigenvalueSolver::computeExpForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallMatrix<double>& expA)
+{
+  EPS eps;
+
+  EPSCreate(PETSC_COMM_SELF, &eps);
+  EPSSetOperators(eps, A, nullptr);
+  EPSSetProblemType(eps, EPS_HEP);
+
+  Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
+
+  PetscInt nb_eigenvalues;
+  EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
+
+  SmallArray<double> eigenvalues(nb_eigenvalues);
+  SmallMatrix<double> P(nb_eigenvalues, nb_eigenvalues);
+  SmallMatrix<double> D(nb_eigenvalues, nb_eigenvalues);
+  D    = zero;
+  expA = SmallMatrix<double>(nb_eigenvalues, nb_eigenvalues);
+  Array<double> eigenvector(nb_eigenvalues);
+  for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
+    Vec Vr;
+    VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr);
+    EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr);
+    VecDestroy(&Vr);
+    for (size_t j = 0; j < eigenvector.size(); ++j) {
+      P(j, i) = eigenvector[j];
+    }
+  }
+
+  for (size_t i = 0; i < eigenvalues.size(); ++i) {
+    D(i, i) = std::exp(eigenvalues[i]);
+  }
+  expA = P * D * transpose(P);
+  EPSDestroy(&eps);
+}
+
 #endif   // PUGS_HAS_SLEPC
 
 EigenvalueSolver::EigenvalueSolver() {}
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
index a9b462509514ace4b736424c65f7c9abe4572f26..febea024a94b7a3b1e033005bf70003a55b9c293 100644
--- a/src/algebra/EigenvalueSolver.hpp
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -23,6 +23,7 @@ class EigenvalueSolver
   void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
                                  SmallArray<double>& eigenvalues,
                                  SmallMatrix<double>& P);
+  void computeExpForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallMatrix<double>& expA);
 #endif   // PUGS_HAS_SLEPC
 
  public:
@@ -63,6 +64,17 @@ class EigenvalueSolver
 #endif   // PUGS_HAS_SLEPC
   }
 
+  template <typename MatrixType>
+  void
+  computeExpForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] SmallMatrix<double>& expA)
+  {
+#ifdef PUGS_HAS_SLEPC
+    this->computeExpForSymmetricMatrix(PETScAijMatrixEmbedder{A}, expA);
+#else    // PUGS_HAS_SLEPC
+    throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
+#endif   // PUGS_HAS_SLEPC
+  }
+
   EigenvalueSolver();
   ~EigenvalueSolver() = default;
 };
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index fb7bf07907bf6c9e1db9333c9713365350918b06..788d1ffd2ba2b7da646820d119a90a27db4272cd 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -114,6 +114,34 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
 #else    //  PUGS_HAS_SLEPC
         REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P),
                             "not implemented yet: SLEPc is required to solve eigenvalue problems");
+#endif   // PUGS_HAS_SLEPC
+      }
+      SECTION("exponential of a matrix")
+      {
+        SmallMatrix<double> expA{};
+        SmallMatrix<double> expS{3};
+
+        expS(0, 0) = 1325.074593930307812;
+        expS(0, 1) = 662.353357244568285;
+        expS(0, 2) = 1324.70671448913637;
+        expS(1, 0) = expS(0, 1);
+        expS(1, 1) = 331.544558063455535;
+        expS(1, 2) = 662.353357244568185;
+        expS(2, 0) = expS(0, 2);
+        expS(2, 1) = expS(1, 2);
+        expS(2, 2) = expS(0, 0);
+
+#ifdef PUGS_HAS_SLEPC
+        EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expA);
+
+        for (size_t i = 0; i < 3; ++i) {
+          for (size_t j = 0; j < 3; ++j) {
+            REQUIRE(expA(i, j) - expS(i, j) == Catch::Approx(0).margin(1E-12));
+          }
+        }
+#else    //  PUGS_HAS_SLEPC
+        REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expP),
+                            "not implemented yet: SLEPc is required to solve eigenvalue problems");
 #endif   // PUGS_HAS_SLEPC
       }
     }