Smoothness of finite-dimensional functional calculus Assume that $f:\mathbb R\to\mathbb R$ is continuous.
Given a real symmetric matrix $A\in\text{Sym}(n)$, we can define $f(A)$ by applying $f$ to its spectrum. More explicitly,
$$ f(A):=\sum f(\lambda)P_\lambda,\qquad A=\sum\lambda P_\lambda. $$
Here both sums are finite, and the second one is the decomposition of $A$ as a linear combination of orthogonal projections ($P_\lambda$ is the projection onto the eigenspace for the eigenvalue $\lambda$, so that $P_\lambda P_{\lambda'}=0$). Such decomposition exists and is unique by the spectral theorem.
I guess it is well known that $f:\text{Sym}(n)\to\text{Sym}(n)$ is continuous.

Assuming $f\in C^\infty(\mathbb R)$, is the induced map $f:\text{Sym}(n)\to\text{Sym}(n)$ also smooth?

I think I can show that it is (Fréchet) differentiable everywhere, but I am wondering whether it is always $C^1$ or even $C^\infty$.
 A: 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.$
A: Yes.  To show that f(A) is n-times differentiable at A=B, simply interpolate f by a polynomial P so that P and its derivatives up to order n agree with f on the spectrum of B.  Clearly P(A) is n-times differentiable, and it isn't too much work to show that f and P have the same nth derivative.
