Analyzing the P.D.E that describes the displacement of a circular membrane The following PDE describes the displacement $u(r,\theta,t)$ of a circular membrane,  $$u_{tt}=c^2\Delta u$$
with B.C. $u(a,\theta,t)=u_r(a,\theta,t)$. I am asked to 
$(a)$ Show that this membrane only oscillates 
$(b)$ Obatain an expression that determines the natural frequencies
$(c)$ Solve the initial value problem if 
$$u(r,\theta,0)=0,\quad u_t(r,\theta,0)=\alpha(r)\sin(3\theta).$$
What I have tried is the following:
After some calculation I was able to seperate the pde into 
$$h''(t)=-\lambda c^2h(t),\quad g''(\theta)=-\mu g(\theta),\quad\text{and}\quad r^2f''(r)+rf'(r)+(\lambda r^2-\mu)f(r)=0$$
where $u(t)=f(r)g(\theta)h(t)$. Also, since the membrane is circular, we have $g(\pi)=g(-\pi)=g'(\pi)=g'(-\pi)$ and $|f(0)|<\infty$.  From the B.C. given we also have $f(a)=-f'(a).$
To answer $(a)$, I assumed we need to show that the solution to $g''=-\mu g$ is a linear combination of sines and cosines, which is trivially shown by using the conditions $g(\pi)=g(-\pi)=g'(\pi)=g'(-\pi)$.
To answer $(b)$, I referred to $h''=-\lambda c^2h$, which with $\lambda>0$, has natural frequencies given by $c\sqrt{\lambda}$.
For $(c)$, using the substitution $z=\sqrt{\lambda}r$ and using the fact that $\mu=n^2$ from the second ode, then we encounter Bessel's equation. Since $|f(0)|<\infty$, then the solution is in the form $f(r)=c_1J_n(\sqrt{\lambda}r)$, however, here I don't know how to implement the fact that $f(a)=f'(a)$. Is my work correct so far? A detailed answer to part $(c)$ would be appreciated. Thank you!
 A: $$
\left|
\begin{array}{l}
u_{tt}=c^2\Delta u\\
u(a,\theta,t)=u_r(a,\theta,t)\\
u(r,\theta,0)=0\\
u_t(r,\theta,0)=\alpha(r)\sin(3\theta)
\end{array}
\right.
$$
I'm only going to focus on the initial/boundary value problem above. Letting some of the eigenvalues be squared is helpful, so I'm going to say that 
$$h''(t)=-\lambda^2 c^2h(t),$$
is the time equation. This lets us write $h(t)=\sin(\lambda c t)$ (up to a scaling constant), which is the solution that vanishes at $t=0$. Similarly, we can see that 
$$
g(\theta) = A\cos(n\theta)+B\sin(n\theta)
$$
meets the the continuity/periodicity conditions for the angular component and solves the equation.
The final piece is the bessel function / radial part, which you noted is solved by 
$$
f(r) = J_n(\lambda r).
$$
To enforce the outward gradient condition on the boundary, we calculate
$$
f(a) = f'(a)
$$
$$
J_n(\lambda a) = \frac{J_{n-1}(\lambda a)-J_{n+1}(\lambda a)}{2}
$$
I claim that for $a$ fixed, for every $n$, the above equation has a family of solutions I'll note as $\lambda_{nm}$ that satisfies this condition. There is no closed form to be had. Here, $m$ can be an integer, and $\lambda_{nm}$ is the $m$th solution of the $n$th equation above. Put another way, $n$ indexes the allowable angular solutions, and $m$ indexes the allowable radial ones, where the radial solutions can only have the discrete $\lambda$ values that satisfy the boundary condition. You can always numerically calculate the exact numbers, but take a look here to see what one example of the roots of this equation looks like:
http://www.wolframalpha.com/input/?i=BesselJ%5B1,x%5D-BesselJ%5B0,x%5D%2F2%2BBesselJ%5B2,x%5D%2F2
Taken together, we get a family of solutions of the form
$$
u_{nm}(r,\theta,t) = A_{nm}\cos(n\theta)\sin(\lambda_{nm} c t)J_n(\lambda_{nm} r)+B_{nm}\sin(n\theta)\sin(\lambda_{nm} c t)J_n(\lambda_{nm} r)
$$
The full general solution is then a sum over the angular and radial modes:
$$
u(r,\theta,t) = \sum_{m=-\infty}^\infty\sum_{n=0}^\infty \left(A_{nm}\cos(n\theta)\sin(\lambda_{nm} c t)J_n(\lambda_{nm} r)+B_{nm}\sin(n\theta)\sin(\lambda_{nm} c t)J_n(\lambda_{nm} r)\right)
$$
To finally solve your IVP, calculate $u_t\left.\right|_{t=0}$
$$
u_t(r,\theta,0) = \sum_{m=-\infty}^\infty\sum_{n=0}^\infty \left(A_{nm}\lambda_{nm} c\cos(n\theta)J_n(\lambda_{nm} r)+B_{nm}\lambda_{nm} c\sin(n\theta)J_n(\lambda_{nm} r)\right)
$$
If this must equal something with only a $\sin(3\theta)$ dependency, that tells us all the $A_{nm}=0$ and only the $n=3$ term in the $n$ summation remains:
$$
u_t(r,\theta,0) = \sin(3\theta)\sum_{m=-\infty}^\infty B_{m}\lambda_{3m} cJ_3(\lambda_{3m} r) = \alpha(r)\sin(3\theta)
$$
Cancelling the $\theta$ dependency
$$
\sum_{m=-\infty}^\infty B_{m}\lambda_{3m} cJ_3(\lambda_{3m} r) = \alpha(r)
$$
For simplicity let $C_m = B_{m}\lambda_{3m} c$ so that
$$
\sum_{m=-\infty}^\infty C_{m} J_3(\lambda_{3m} r) = \alpha(r)
$$
Now the task is to find the $C_{m}$ coefficients that make this true, i.e., find a representation of an arbitrary $\alpha$ function in terms of $J_3$ bessel functions with the set of discrete frequencies determined by the radial boundary conditions. Again, this can all be done numerically to any precision desired, but there is no closed form to be had here. Once you find those coefficients, the final solution will be
$$
u(r,\theta,t) = \sin(3\theta)\sum_{m=-\infty}^\infty B_{m}\sin(\lambda_{3m} c t)J_3(\lambda_{3m} r)
$$
