No, an "analytical method (without solve a system of equations)" for solving such equations does not exist. This is because any matrix equation is by definition a system of equations.
However, considering when a matrix is already diagonalized, finding the solution is very simple. Read ahead for some considerations of the non-diagonalizable case.
Every matrix may be orthogonaly reduced to triangular. Take then the equation $Y^2 = R$ where $R$ is triangular and $C=QRQ^*$ where $QQ^*=I$. Then the solution for $X$ may be found
$$X^2=C=QRQ^* = QY^2Q^*$$
So that
$$ X = QYQ^*$$
So it comes down to finding for a triangular matrix. Lets look at the $3 \times 3$ case.
$$\pmatrix{d_0 & 0 & 0 \\ a_1 & d_1 & 0 \\ a_2 & b_2 & d_2}\pmatrix{d_0 & 0 & 0 \\ a_1 & d_1 & 0 \\ a_2 & b_2 & d_2} = \pmatrix{d_0^2 & 0 & 0 \\ a_1(d_0 + d_1) & d_1^2 & 0 \\ a_2(d_0 + d_2) + a_1b_2 & b_2(d_1 + d_2) & d_2^2}$$
This would not necessarily give a result " analytical method (without solve a system of equations)" but it does show that the square root does exist, and there are $2^n$ different solutions, where $n$ is the number of non-zero eigenvalues.
EDIT: I will rollback if this is wrong as it is ripped from a text, but...
Functions of Jordan Blocks
For a $k \times k$ Jordan block $J_*$ with eigenvalue $\lambda$, and for a function $f(z)$ such that $f(\lambda), f'(\lambda), \dots, f^{(k-1)}(\lambda)$ exist, $f(J_*)$ is defined to be
$$f(J_*) = f\pmatrix{\lambda & 1 \\ & \ddots & \ddots \\ & & \ddots & 1 \\ &&& \lambda}=\pmatrix{f(\lambda) & f'(\lambda) & \frac{f''(\lambda)}{2!} & \cdots & \frac{f^{(k-1)}(\lambda)}{(k-1)!}\\ & f(\lambda)& f'(\lambda)& \ddots & \vdots\\ & & \ddots & \ddots & \frac{f''(\lambda)}{2!}\\ & & & f(\lambda) & f'(\lambda)\\ & & & & f(\lambda)\\}$$