diff --git a/code/math.f90 b/code/math.f90
index 7e090de02..47c7da175 100644
--- a/code/math.f90
+++ b/code/math.f90
@@ -1936,8 +1936,8 @@ end subroutine math_spectralDecompositionSym
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief eigenvalues and eigenvectors of symmetric 3x3 matrix m using an analytical expression
-!> and the general LAPACK powered version as fallback
+!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m using an analytical expression
+!> and the general LAPACK powered version for arbritrary sized matrices as fallback
 !> @author Joachim Kopp, Max–Planck–Institut für Kernphysik, Heidelberg (Copyright (C) 2006)
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
 !> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3)
@@ -1958,7 +1958,7 @@ subroutine math_spectralDecompositionSym33(m,values,vectors)
                     m(1, 2)**2_pInt]
 
  T = maxval(abs(values))
- U = MAX(T, T**2_pInt)
+ U = max(T, T**2_pInt)
  threshold = sqrt(5.0e-14_pReal * U**2_pInt)
 
 ! Calculate first eigenvector by the formula v[0] = (m - lambda[0]).e1 x (m - lambda[0]).e2
@@ -2051,35 +2051,35 @@ end function math_eigenvaluesSym
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Eigenvalues of symmetric 33 matrix m
+!> @brief eigenvalues of symmetric 33 matrix m using an analytical expression
+!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @details similar to http://arxiv.org/abs/physics/0610206 (DSYEVC3)
+!> but apparently more stable solution and has general LAPACK powered version for arbritrary sized
+!> matrices as fallback
 !--------------------------------------------------------------------------------------------------
 function math_eigenvaluesSym33(m)
 
  implicit none
  real(pReal), intent(in), dimension(3,3) :: m
- real(pReal), dimension(3) :: math_eigenvaluesSym33, invariants
- real(pReal) :: R, S, T,  P, Q, rho, phi
+ real(pReal), dimension(3) :: math_eigenvaluesSym33,invariants
+ real(pReal) :: P, Q, rho, phi
  real(pReal), parameter :: TOL=1.e-14_pReal
 
  invariants = math_invariantsSym33(m)
 
- R=-invariants(1)
- S= invariants(2)
- T=-invariants(3)
+ P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal
+ Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3)
 
- P=S-R**2.0_pReal/3.0_pReal
- Q=2.0_pReal/27.0_pReal*R**3.0_pReal-R*S/3.0_pReal+T
-
- if((abs(P) < TOL) .and. (abs(Q) < TOL)) then
-   math_eigenvaluesSym33 = invariants(1)/3.0_pReal
+ if(any(abs([p,q]) < TOL)) then
+   math_eigenvaluesSym33 = math_eigenvaluesSym(m)
  else
    rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
-   phi=acos(math_limit(-Q/rho/2.0_pReal,-1.0_pReal,1.0_pReal))
+   phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
    math_eigenvaluesSym33 = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
                          [cos(phi/3.0_pReal), &
-                          cos(phi/3.0_pReal+2.0_pReal/3.0_pReal*PI), &
-                          cos(phi/3.0_pReal+4.0_pReal/3.0_pReal*PI) &
-                         ] -R/3.0_pReal
+                          cos((phi+2.0_pReal*PI)/3.0_pReal), &
+                          cos((phi+4.0_pReal*PI)/3.0_pReal) &
+                         ] + invariants(1)/3.0_pReal
  endif
 end function math_eigenvaluesSym33