Smoothness of finite-dimensional functional calculus

Yes. The can be derived from the resolvent formalism.

I'll just do the $C^1$ case and leave higher derivatives as an exercise - ask if it's not clear how to generalize. I am basically using formula (2.7) of "On differentiability of symmetric matrix valued functions", Alexander Shapiro, http://www.optimization-online.org/DB_HTML/2002/07/499.html. (I think there's a typo in the middle case of the display preceeding (2.7): $f(\mu_j)/(\mu_j-\mu_k)$ should be $(f(\mu_j)-f(\mu_k))/(\mu_j-\mu_k).$) Shapiro's paper references "(cf., [4])", where [4] is the 600-page textbook "Perturbation Theory for Linear Operators" by T. Kato, but I don't know if that is helpful for this specific question.

I will call the induced map $f^*$ to distinguish it from $f.$ I'll also call the dimension $p$ instead of $n.$

It suffices to show $f^*$ is $C^1$ for matrices with eigenvalues in a given bounded interval $J.$ Approximate $f$ by polynomials $f_n$ such that $\sup_{x\in J}|f(x)-f_n(x)|\to 0$ and $\sup_{x\in J}|f'(x)-f_n'(x)|\to 0.$ Since $f_n$ is analytic, $f^*_n$ can be evaluated using resolvents:

$$f_n^*(X) = \frac{1}{2\pi i}\int_C f_n(z)(z I_p - X)^{-1} dz$$ where $C$ is an anticlockwise circle in the complex plane with $J$ in its interior. For $H\in\mathrm{Sym}(p),$

\begin{align*} f_n^*(X+H) &= \frac{1}{2\pi i}\int_C f_n(z)(z I_p - X-H)^{-1} dz\\ &= \frac{1}{2\pi i}\int_C f_n(z)(z I_p - X)^{-1}+f_n(z)(z I_p - X)^{-1}H(z I_p - X)^{-1} +\dots dz\\ &= \frac{1}{2\pi i}\int_C f_n(z)\sum_{\lambda}(z-\lambda)^{-1}P_\lambda +f_n(z)\sum_{\lambda_1,\lambda_2}(z-\lambda_1)^{-1}(z-\lambda_2)^{-1}P_{\lambda_1}HP_{\lambda_2}+\dots dz\\ &= f_n^*(X)+\sum_{\lambda_1,\lambda_2} P_{\lambda_1} H P_{\lambda_2}\int_0^1 f'_n(t\lambda_1+(1-t)\lambda_2)+\dots dt \end{align*}

The second equality uses the Taylor expansion $$(A-H)^{-1}=A^{-1}+A^{-1}HA^{-1}+\dots$$ with $A=z I_p-X.$ The third equality uses $(zI_p - X)^{-1}=\sum_\lambda (z-\lambda)^{-1} P_\lambda.$ The fourth equality uses $\int_C f_n(z)(z-\lambda)^{-1}(z-\mu)^{-1}dz =\int_0^1 f'_n(t\lambda+(1-t)\mu)dt.$

This gives a bound

$$\|Df^*_n(X)H\| \leq c_p\|H\|\cdot \sup_{x\in J}|f'_n(x)-f_n'(x)|$$

for some constant $c_p>0,$ where $\|\cdot\|$ is any matrix norm. This shows that $f^*$ can be approximated arbitrarily well in the $C^1$ norm, which means it's $C^1.$