Following Pete's advice, here is my comment with some more details added:

For a Jordan block $A$ with eigenvalue $t\neq 0$ write $A=tB$. $B$ has 1's on the main diagonal and $1/t$'s immediately above it. $N=B-I$ is nilpotent, so a square root $C$ of $B=I+N$ can be found using the binomial formula (which gives a finite sum). Then $\sqrt{t}C$ will be a square root of $A$.

If $A$ is arbitrary, then find the Jordan form $B$ of $A$ so that $A=C^{-1}BA$. If there are no zero eigenvalues, then we can find a square root of each block and then conjugate back.

If there are Jordan blocks with eigenvalue 0, the problem gets a bit trickier. The square of a Jordan $m$ by $m$ block with zero eigenvalue is conjugate to the union of two $m/2$ by $m/2$ blocks if $m$ is even and to the union of an $(m-1)/2$ by $(m-1)/2$ block and an $(m+1)/2$ by $(m+1)/2$ block if $m$ is odd. This allows one to compute a square root of a union of two Jordan blocks of equal sizes or of a union of an $n$ by $n$ block and an $(n+1)$ by $(n+1)$ block.

Let $a_1\leq \ldots\leq a_k$ be the sizes of the zero eigenvalue Jordan blocks (including 1 by 1 ones, so $\sum a_i$ is the dimension of the generalized eigenspace with eigenvalue 0). $A$ has a square root, iff $a_1\ldots,a_k$ can be obtained from a sequence $b_1\leq \ldots\leq b_l$ of positive integers by replacing an even $m$ with $m/2,m/2$, an odd $m$ with $(m-1)/2,(m+1)/2$ and leaving 1's untouched. (I know this looks messy but can't think of anything better.)

Of course, a square root of a matrix is not unique (if it exists).

Note that if all eigenvalues of $A$ are positive, then numerically it's probably easier to use the binomial formula straight away:

$$\sqrt{A}=\sqrt{t}(I+\frac{1}{2t}X-\frac{1}{8t}X^2 +\cdots).$$

Here $X=A-tI$ and the formula is valid for $t$ greater then the maximum eigenvalue of $A$.

Following Pete's advice, here is my comment with some more details added:

For a Jordan block $A$ with eigenvalue $t\neq 0$ write $A=tB$. $B$ has 1's on the main diagonal and $1/t$'s immediately above it. $N=B-I$ is nilpotent, so a square root $C$ of $B=I+N$ can be found using the binomial formula (which gives a finite sum). Then $\sqrt{t}C$ will be a square root of $A$.

If $A$ is arbitrary, then find the Jordan form $B$ of $A$ so that $A=C^{-1}BA$. If there are no zero eigenvalues, then we can find a square root of each block and then conjugate back.

If there are Jordan blocks with eigenvalue 0, the problem gets a bit trickier. The square of a Jordan $m$ by $m$ block with zero eigenvalue is conjugate to the union of two $m/2$ by $m/2$ blocks if $m$ is even and to the union of an $(m-1)/2$ by $(m-1)/2$ block and an $(m+1)/2$ by $(m+1)/2$ block if $m$ is odd. This allows one to compute a square root of a union of two Jordan blocks of equal sizes or of a union of an $n$ by $n$ block and an $(n+1)$ by $(n+1)$ block.

Let $a_1\leq \ldots\leq a_k$ be the sizes of the zero eigenvalue Jordan blocks (including 1 by 1 ones, so $\sum a_i$ is the dimension of the generalized eigenspace with eigenvalue 0). $A$ has a square root, iff $a_1\ldots,a_k$ can be obtained from a sequence $b_1\leq \ldots\leq b_l$ of positive integers by replacing an even $m$ with $m/2,m/2$, an odd $m$ with $(m-1)/2,(m+1)/2$ and leaving 1's untouched. (I know this looks messy but can't think of anything better.)

Of course, a square root of a matrix is not unique (if it exists).

Note that if all eigenvalues of $A$ are positive, then numerically it's probably easier to use the binomial formula straight away:

$$\sqrt{A}=\sqrt{t}(I+\frac{1}{2t}X-\frac{1}{8t^2}X^2 +\cdots).$$

Here $X=A-tI$ and the formula is valid for $t$ greater then the maximum eigenvalue of $A$.