I'm trying to calculate the value of a matrix function. As far as I understood, this is done by first decomposing my matrix $A$ into $PJP^{-1}$. Where $J$ is in Jordan normal form.
However, this appears to be a ton of work for $3\times3$ matrices so I'm thinking I might be using the wrong algorithm.
The way I do it is first calculating the characteristic polynomial. This is already one computation of a $3\times3$ determinant and then simplification. After than, I solve the $(A-\lambda I)x=0$ for each eigenvalue to get the geometric multiplicities. Then, as it usually is the case, I end up missing 2 of 3 eigenvectors so I don't have a Jordan basis. This means that I need to solve the $AP=PJ$ system where I usually only know the first column of $P$. After all this is done, I have to find the inverse of $P$.
Then, solving the task at hand is relatively simple and requires computing a couple of derivatives and a few matrix multiplications.
How can this process be simplified and sped up? What are the techniques I can apply here in general? Are there any that can be applied for some specific cases of $A$?