A vibrating drumhead can be modeled by wave equation in polar coordinates and appropriate boundary and initial conditions. The modeling problem can be set up as follows. We assume that our drumhead is a round disk with unit radius ($r=1$).

PDE: $u_{tt}=c^2\left(u_{rr}+\frac{1}{r}u_r+\frac{1}{r^2}u_{\theta\theta}\right),\ 0<r<1$

BC: $u=0$ when $r=1$, $0<t<\infty$

ICs: \begin{align*}u(r,\theta,0)&=f(r,\theta),\\u_t(r,\theta,0)&=g(r,\theta)
\end{align*}
Note that the expression inside () is the Laplacian $\nabla^2 u$ in polar coordinates, i.e. $$\nabla^2u=u_{rr}+\frac{1}{r}u_r+\frac{1}{r^2}u_{\theta\theta}.$$
Let $u(r,\theta,t)=U(r,\theta)T(t)$. $U(r,\theta)$ is the shape of the drumhead while $T(t)$ is the oscillatory factor. Then the PDE can be broken into two differential equtions, one PDE and one ODE:
\begin{align*}
\nabla^2U+\lambda^2U&=0\ (\mbox{Helmholtz equation})\\\ddot{T}+\lambda^2c^2T&=0\ (\mbox{Simple harmonic motion})
\end{align*}
With the boundary condition, we have the Helmholtz eigenvalue problem:
\begin{align*}
\nabla^2U+\lambda^2U&=0\\
U(1,\theta)&=0.
\end{align*}
To solve this Helmholtz eigenvalue problem, we assume that $U(r,\theta)=R(r)\Theta(\theta)$. Then the Helmholtz equation is broken into two ODEs with BCs:
\begin{align*}
r^2R^{\prime\prime}+rR’+(\lambda^2r^2-n^2)R&=0\ (\mbox{Bessel’s equation}),\\
R(1)&=0
\end{align*}
and
\begin{align*}
\Theta^{\prime\prime}+n^2\Theta&=0,\\
\Theta(0)&=\Theta(2\pi).
\end{align*}
The boundary condition for $\Theta(\theta)$ is imposed because we want $\Theta$ to be periodic with period $2\pi$. We also impose an additional condition for $R(r)$:
$$R(0)<\infty\ (\mbox{physical condition}).$$

Let us discuss the solutions of the Bessel’s equation. The Bessel’s equation has two independent solutions
$J_n(\lambda r)$ and $Y_n(\lambda r)$, where $$J_n(x)=\sum_{k=0}^\infty\frac{(-1)^kx^{2k+n}}{2^{2k+n}k!(n+k)!}$$ and $$Y_n(x)=\frac{2}{\pi}J_n(x)\left[\ln\left(\frac{x}{2}\right)+\gamma\right]+\frac{2^n}{\pi x^n}\sum_{k=0}^\infty\frac{\beta_{mk}}{2^{2k}k!}x^{2k}.$$
Here, $\gamma$ is called the Euler-Macheroni constant and is given by
$$\gamma=\lim_{n\to\infty}\left[\left(\sum_{k=1}^n\frac{1}{k}\right)-\ln n\right]$$
and $\beta_{mk}$ is defined by
$$\beta_{mk}=\left\{\begin{array}{ccc}
-(n-1-k)! & \rm{if} & k\leq n-1,\\(-1)^{k-n-1}\frac{h_{k-n}+h_k}{(k-n)!} & \rm{if} & k>m-1,
\end{array}\right.$$ where $h_p=\sum_{i=1}^p\frac{1}{i}$.
$J_n(x)$ is called Bessel function of the first kind and $Y_n(x)$ is Bessel function of the second kind. $Y_n(x)$ is also called Neumann function. The general solution $R(r)$ of the Bessel’s equation is then written as
$$R(r)=AJ_n(\lambda r)+BY_n(\lambda r).$$ Since $Y_n(\lambda r)$ is not defined at $r=0$, due to the boundary condition $R(0)<\infty$, we require that $B=0$. So we consider only Bessel function of the first kind $J_n(x)$ here. Let us plot Bessel function of the first kind using Maxima. In Maxima, Bessel function of the first kind of order $n$ and argument $x$, $J_n(x)$ is denoted by bessel_j(n,x). To plot $J_1(x)$, $0\leq x\leq 10$, type the command

plot2d(bessel_j(1,x),[x,0,10]);

and press enter as you see in the following screenshot.

Then you will see the plot.

If you want to plot $J_2(x)$, type the command

plot2d(bessel_j(2,x),[x,0,10]);

and press enter.

If you want to plot multiple Bessel functions, for example,  $J_1(x)$, $J_2(x)$, $J_3(x)$, $J_4(x)$,  type the command

plot2d([bessel_j(1,x),bessel_j(2,x),bessel_j(3,x),bessel_j(4,x)],[x,0,10]);

and press enter.