There is also an unconditional deterministic polynomial-time algorithm to find $x,y,z,w \in \mathbf{F}_p^\times$ such that $x^2+y^2+z^2+w^2=n$, given any $n \in \mathbf{Z}$ and any prime $p \ge 7$.
First, given $a,b,c \in \mathbf{F}_p^\times$, Theorem 1.10 of Christiaan van de Woestijne's thesis lets one find an $\mathbf{F}_p$-point $P$ on the smooth conic $C \colon ax^2+by^2=cz^2$ in $\mathbf{P}^2$ over $\mathbf{F}_p$. The usual trick of drawing lines through $P$ and taking the second point of intersection with $C$ lets one parametrize $C(\mathbf{F}_p)$. At most $6$ points of $C(\mathbf{F}_p)$ have one of $x,y,z$ equal to $0$, so by trying at most $7$ lines, one finds a point on the affine curve $ax^2+by^2=c$ with $x,y \in \mathbf{F}_p^\times$.
Now, to solve $x^2+y^2+z^2+w^2=n$, choose $c \in \mathbf{F}_p \setminus \{0,n\}$, and apply the previous sentence to find $x,y,z,w \in \mathbf{F}_p^\times$ satisfying $x^2+y^2=c$ and $z^2+w^2=n-c$.
Here's a copy-paste of something I wrote up a while ago:
Lemma: Let $q$ be odd, and let $Q$ be the set of quadratic residues (including $0$) in $\mathbb F_q$. Then the number of elements $s_q(c)$ in $\{x^2+c|x \in \mathbb{F}_q\} \cap Q$ is given by
\begin{array}{|c|c|c|}
\hline
& c \in Q & c \notin Q \\
\hline
-1 \in Q & \frac{q+3}{4} & \frac{q-1}{4} \\
\hline
-1 \notin Q & \frac{q+1}{4} & \frac{q+1}{4} \\
\hline
\end{array}
Proof: If, for $x,y,c\in \mathbb{F}_q,\ c \neq 0$ we have $x^2+c=y^2$, then $c=y^2-x^2=(y-x)(y+x)$. Now for all the $q-1$ elements $d\in \mathbb{F}_q^{\ast}$, we can let $y-x=d$ and $y+x=\frac{c}{d}$. But the pairs $(d,\frac{c}{d}),(-d,\frac{c}{-d}),(\frac{c}{d},d),(\frac{c}{-d},-d)$ all give the same value of $y^2=\frac{1}{4}(d+c/d)^2$. Also, as $q$ is odd, $d\neq -d\ \forall d$. But if $c\in Q$, for $2$ values of $d$ we have $d=\frac{c}{d}$ and if $-c\in Q$, for 2 values of $d$ we have $d=\frac{c}{-d}$. So we have
$$ s_q(c) = \left\{ \begin{array}{rcll}
\frac{\frac{q-1}{2}-2}{2}+2 & = & \frac{q+3}{4} & if\ c\in Q,\ -c\in Q \\
\frac{\frac{q-1}{2}-1}{2}+1 & = & \frac{q+1}{4} & if\ c\in Q,\ -c\notin Q \\
\frac{\frac{q-1}{2}-1}{2}+1 & = & \frac{q+1}{4} & if\ c\notin Q,\ -c\in Q \\
& & \frac{q-1}{4} & if\ c\notin Q,\ -c\notin Q
\end{array} \right. $$
and hence the result.
Best Answer
This probability can be calculated exactly, and indeed it approaches $1/2$ rather quickly — more precisely, for each $p$ it approaches the fraction $(p-1)/(2p)$ of quadratic residues $\bmod p$. This can be proved by elementary means, but perhaps the nicest way to think about it is that if you choose $n$ numbers $a_i$ independently and sum $a_i^2 \bmod p$, the resulting distribution is the $n$-th convolution power of the distribution of a random single square — so its discrete Fourier transform is the $n$-th power of the D.F.T., call it $\gamma$, of the distribution of $a^2 \bmod p$. For this purpose $\gamma$ is normalized so $\gamma(0)=1$. Then for $k \neq 0$ we have $\gamma(k) = (k/p) \gamma(1)$ [where $(\cdot/p)$ is the Legendre symbol], and $$ p \gamma(1) = \sum_{a \bmod p} \exp(2\pi i a^2/p), $$ which is a Gauss sum and is thus a square root of $\pm p$. It follows that $|\gamma(k)| = p^{-1/2}$, from which we soon see that each value of the convolution approaches $1/p$ at the exponential rate $p^{-n/2}$, and the probability you asked for approaches $(p-1)/(2p)$ at the same rate.
As noted above, this result, and indeed the exact probability, can be obtained by elementary means, yielding a (known but not well-known) alternative proof of Quadratic Reciprocity(!). But that's probably too far afield for the present purpose.