Rrerku

Assume that $f:mathbb Rtomathbb R$ is continuous.
Given a real symmetric matrix $AintextSym(n)$, we can define $f(A)$ by applying $f$ to its spectrum. More explicitly,
$$ f(A):=sum f(lambda)P_lambda,qquad A=sumlambda 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:textSym(n)totextSym(n)$ is continuous.

Assuming $fin C^infty(mathbb R)$, is the induced map $f:textSym(n)totextSym(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$.

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_xin J|f(x)-f_n(x)|to 0$ and $sup_xin J|f'(x)-f_n'(x)|to 0.$ Since $f_n$ is analytic, $f^*_n$ can be evaluated using resolvents:

$$f_n^*(X) = frac12pi iint_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 $HinmathrmSym(p),$

beginalign*
f_n^*(X+H)
&= frac12pi iint_C f_n(z)(z I_p - X-H)^-1 dz\
&= frac12pi iint_C f_n(z)(z I_p - X)^-1+f_n(z)(z I_p - X)^-1H(z I_p - X)^-1 +dots dz\
&= frac12pi iint_C f_n(z)sum_lambda(z-lambda)^-1P_lambda +f_n(z)sum_lambda_1,lambda_2(z-lambda_1)^-1(z-lambda_2)^-1P_lambda_1HP_lambda_2+dots dz\
&= f_n^*(X)+sum_lambda_1,lambda_2 P_lambda_1 H P_lambda_2int_0^1 f'_n(tlambda_1+(1-t)lambda_2)+dots dt
endalign*

The second equality uses the Taylor expansion $$(A-H)^-1=A^-1+A^-1HA^-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)^-1dz =int_0^1 f'_n(tlambda+(1-t)mu)dt.$

This gives a bound

$$|Df^*_n(X)H| leq c_p|H|cdot sup_xin 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.$

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.