We're looking for solutions $u(\mathbf{r},t)$ to the wave equation
$$\frac{\partial^2 u}{\partial t^2}=c^2\nabla^2u$$
Which, in spherical coordinates, is
$$\frac{c^2}{r^2\sin(\phi)}\left(\frac{\partial }{\partial r}\left(r^2\sin(\phi)\frac{\partial u}{\partial r}\right)+\frac{\partial }{\partial \theta}\left(\frac{1}{\sin(\phi)}\frac{\partial u}{\partial \theta}\right)+\frac{\partial }{\partial \phi}\left(\sin(\phi)\frac{\partial u}{\partial \phi}\right)\right)=\frac{\partial^2u}{\partial t^2}$$
Using the $(r,\theta,\phi)$ (radius, azimuthal angle, polar angle) convention. First step: Assume the solution is separable, i.e, $u(\mathbf{r},t)=R(r)\Theta(\theta)\Phi(\phi)T(t)$. Now we plug this into our wave equation.
$$R\Theta\Phi T''=\frac{c^2}{r^2\sin\phi}\left(\frac{\partial }{\partial r}\left(r^2\sin(\phi) R'\Theta\Phi T\right)+\frac{\partial }{\partial \theta}\left(\frac{1}{\sin\phi}R\Theta'\Phi T\right)+\frac{\partial }{\partial \phi}\left(\sin(\phi) R\Theta\Phi'T\right)\right)$$
Simplifying,
$$R\Theta\Phi T''=c^2\left(\frac{\Theta\Phi T}{r^2}\frac{\partial }{\partial r}\left(r^2R'\right)+\frac{R\Phi T}{r^2\sin^2\phi}\frac{\partial }{\partial \theta}\left(\Theta'\right)+\frac{R\Theta T}{r^2\sin\phi}\frac{\partial }{\partial \phi}\left(\sin(\phi) \Phi'\right)\right)$$
Dividing out,
$$ \frac{T''}{T}=\frac{c^2}{r^2}\left(\frac{1 }{R}\frac{\partial }{\partial r}\left(r^2R'\right)+\frac{1 }{\Theta \sin^2\phi}\frac{\partial }{\partial \theta}\left(\Theta'\right)+\frac{1 }{\Phi \sin\phi}\frac{\partial }{\partial \phi}\left(\sin(\phi) \Phi'\right)\right)$$
Now I've illustrated all this so far since this is how you'd approach this problem in general. However, since the initial conditions you've provided are radially symmetric, i.e, independent of $\theta$ and $\phi$, we can assume $\Theta'=\Phi'=0.$ Therefore the above simplifies to
$$\frac{T''}{T}=\frac{c^2}{r^2}\frac{1 }{R}\frac{\partial }{\partial r}\left(r^2R'\right)$$
We now assume our solution is "pure frequency", i.e,
$$\frac{T''}{T}=-\omega^2 \implies T(t)=a_1\cos(\omega t)+a_2\sin(\omega t)$$
Substituting into the previous,
$$-\omega^2=\frac{c^2}{r^2}\frac{1 }{R}\frac{\partial }{\partial r}\left(r^2R'\right)$$
Defining $k=\frac{\omega}{c}$,
$$\frac{\partial }{\partial r}\left(r^2R'\right)+k^2r^2R=0$$
This is an Emden-Fowler differential equation, and has the solution
$$R(r)=\frac{1}{r}\left(b_1e^{-ikr}+\frac{ib_2}{2k}e^{ikr}\right)$$
Let $b_1=p+iq$, $b_2=p'+iq'$. Then,
$$r\cdot R(r)=(p+iq)e^{-ikr}+\frac{i(p'+iq')}{2k}e^{ikr}$$
$$r\cdot R(r)=(p+iq)e^{-ikr}+\frac{1}{2k}(-q'+ip')e^{ikr}$$
Using Euler's formula,
$$r\cdot R(r)=p\cos(kr)-ip\sin(kr)+iq\cos(kr)+q\sin(kr)+\frac{1}{2k}(-q'\cos(kr)+iq'\sin(kr)+ip'\cos(kr)-p'\sin(kr))$$
Grouping the cosine and sine terms together,
$$r\cdot R(r)=\cos(kr)\left(p+iq-\frac{q'}{2k}+\frac{ip'}{2k}\right)+\sin(kr)\left(-ip+q+\frac{iq'}{2k}-\frac{p'}{2k}\right)$$
In order for $R(r)$ to be real-valued, we require
$$q+\frac{p'}{2k}=0 ~\text{ and } -p+\frac{q'}{2k}=0$$
Thus,
$$\begin{bmatrix}
p'\\
q'
\end{bmatrix} =2k\begin{bmatrix}
0 & -1\\
1 & 0
\end{bmatrix}\begin{bmatrix}
p\\
q
\end{bmatrix}$$
Therefore given $b_1=p+iq$, we require $b_2=2k(-q+ip)$. Thus,
$$R(r)=\frac{1}{r}\left(\cos(kr)(p-\frac{2kp}{2k})+\sin(kr)(q-\frac{-2kq}{2k})\right)$$
Finally, substituting $k=\frac{\omega}{c}$ back in and simplifying,
$$R(r)=2q\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)$$
With $\operatorname{sinc}(x):=(1/x)\sin(x)$. Side note: $\operatorname{sinc}$ is the zeroth spherical Bessel function of the first kind, knowledge of which will be needed both for the more general PDE involving $\theta$ and $\phi$, and dealing with the initial conditions later on.
Now, our solution $u(\mathbf{r},t)$(which we can restate as simply $u(r,t)$ because of no $\theta$ or $\phi$ dependence) for this radially symmetric wave equation is
$$u(r,t)=\left(a_1\cos(\omega t)+a_2\sin(\omega t)\right)2q\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)$$
Or, more nicely,
$$u(r,t)=C\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)(A\cos(\omega t)+B\sin(\omega t))$$
Now let's examine our initial conditions. We want $\frac{\partial u}{\partial t}(r,0)=0$.
$$\frac{\partial u}{\partial t}=C\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)(-A\sin(\omega t)+B\cos(\omega t))$$
$$\frac{\partial u}{\partial t}\bigg|_{(r,0)}=C\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)B$$
In order for this to be non trivially $=0 ~ \forall r$, we need $B=0$. As both $A$ and $C$ are constants we can rename $A\cdot C \to A$, therefore
$$u(r,t)=A\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)\cos(\omega t)$$
Before we deal with the other initial condition, we need to generalize our solution a bit.
Recall that we can state $u(r,t)$ in terms of spherical Bessel functions:
$$u(r,t)=A\frac{\omega}{c}\cos(\omega t)j_0\left(\frac{\omega}{c}r\right)$$
Or, better, in terms of ordinary Bessel functions using the identities here:
$$u(r,t)=A\frac{\omega}{c}\cos(\omega t)\sqrt{\frac{\pi}{2r}}J_{1/2}\left(\frac{\omega}{c}r\right)$$
Now, one thing we can notice is that the wave equation is linear, i.e if $u_1$ and $u_2$ are both solutions, $a\cdot u_1 +b\cdot u_2$ is also a solution, given real numbers $a$ and $b$. Therefore we can restate our solution more generally as a linear combination of pure frequency solutions (also known as normal modes),
$$u(r,t)=\frac{A}{c}\sqrt{\frac{\pi}{2r}}\sum_{n=0}^\infty a_n\cos(\omega_nt)\omega_n J_{1/2}\left(\frac{\omega_n}{c}r\right)$$
Where $a_1,a_2,...$ is any sequence of real numbers (provided the sum converges), and $\omega_n$ is a factor that ensures we have roots at $r=\pm 1$. Specifically, let $\alpha_{1/2,n}$ be the $n$th positive root of $J_{1/2}$. Then $\omega_n =c \alpha_{1/2,n}$. Also let $A_n=A a_n\alpha_n$. We can now state $u(r,0)$ as
$$u(r,0)=\sqrt{\frac{\pi}{2r}}\sum_{n=0}^\infty A_n J_{1/2}\left(\alpha_{1/2,n} r\right)$$
Recall that in this example, our initial conditions are
$$ u(r,0)=u_0(r) = \begin{cases} \sqrt{1-r^2}, & \text{if } |r| \leq 1, \\ 0 & \text{if } |r| \geq 1 \end{cases} $$
So now the challenge is to find a sequence of real numbers $A_1,A_2,...$ such that
$$u_0(r)\sqrt{\frac{2r}{\pi}}=\sum_{n=0}^\infty A_n J_{1/2}\left(\alpha_{1/2,n} r\right)$$
For which you'll need Fourier-Bessel series and/or Hankel Transforms. In this case it turns out that
$$A_n=\frac{2}{J_{3/2}(\alpha_{1/2,n})^2}\int_0^1 r\sqrt{\frac{2r}{\pi}}\sqrt{1-r^2} J_{1/2}(\alpha_{1/2,n}r)\mathrm{d}r$$
If you can't be bothered with that, you can use the very rough approximation
$$\sqrt{1-x^2}\approx \frac{\sin(2x)}{2x}=j_0(2x)=\sqrt{\frac{\pi}{2x}}J_{1/2}(2x).$$
I know this is not a complete answer, but I hope it's a good start!
