I want to solve the following system of differential equations $$\dot x_0(t)=-ax_0(t)+\mu x_1(t)$$ $$\dot x_n(t)=(\lambda(n+1)+a)x_{n-1}(t)-(a+(\lambda+\mu)n)x_n(t)+\mu(n+1)x_{n+1} \qquad \forall n\geq 1.$$
To do this I wanted to calculate the matrix exponential $\exp(tA)$ where A is the matrix associated with the given system. Unfortunately I wasn't able to do this. Is this even the right strategy?
I don't know the origin of your set of equations, but an infinite dimensional dynamical system smells a lot like a PDE to me. As it turns out, you can convert the above to a (linear!) PDE as follows. Introduce \begin{equation} y(z,t) = \sum_{n=0}^\infty x_n(t)\,z^n. \end{equation} The recursive relation you posted is equivalent to demanding that $y$ satisfies \begin{equation} \frac{\partial y}{\partial t} = \left(z(\lambda+a)-a\right)\,y + (1-\lambda-\mu)z \frac{\partial y}{\partial z} + \mu \frac{\partial y}{\partial z}, \end{equation} if I haven't made any mistakes. This linear first order PDE can now be solved using the method of characteristics. The series expansion of the solution in powers of $z$ will give you $x_n(t)$; to find a specific $x_n$, calculate $\frac{1}{n!}\frac{\partial^n y}{\partial z^n}(0,t)$. FYI: This could be called the 'generating function approach'.