How can I solve the wave equation for a circular rod in cylindrical coordinates using finite differences?

Question

I have a problem with the stability of finite difference method for the wave equation in cylindrical coordinates.

the equation is: $$ \frac{\partial^2 \omega_n}{\partial r^2}+\frac{1}{r}\frac{\partial \omega_n}{\partial r}+\frac{\partial^2\omega_n}{\partial z^2}-\frac{\omega_n}{r^2}=\frac{1}{c^2}\frac{\partial^2\omega_n}{\partial t^2} $$

and when I run it in Matlab, I get the following result:

Does anyone know where the problem is?

I put the code here so that if anyone knows where the problem is, help me out.

%domain
clear;
Lx=10;
Ly=10;
dx=0.1;
dy=dx;
nx=fix(Lx/dx);
ny=fix(Ly/dy);
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);
%time
T=10;
%varables
wn=zeros(nx,ny);
wnm1=wn; %w at time n-1
wnp1=wn; %w at tome n+1
%parameters
CFL=0.5; %CFL=c.dt/dx
c=1;
dt=(CFL*dx/c);

%dt=1e-6;

%%initial conditions
%%time stepping loop
t=0;
while(t<T)
    %reflecting boundary conditions
    wn(:,[1 end])=0;
    wn([1 end],:)=0;


    %solution
    t=t+dt;
    wnm1=wn; wn=wnp1; %save current and previous arrays

    %sourse
    wn(50,50)=dt^2*20*sin(30*pi*t/20);
    for i=2:nx-1 for j=2:ny-1

           wnp1(i,j)=2*wn(i,j)-wnm1(i,j)...
               +CFL^2*(((wn(i+1,j)-2*wn(i,j)+wn(i-1,j))/dx^2+(1/(2*i*dx^2))*(wn(i+1,j)-wn(i-1,j))+(wn(i,j+1)-2*wn(i,j)+wn(i,j-1))/dy^2-wn(i,j)/(i*dx)^2));

        end,end
    %chek convergence

    %visualize at selected steps
    clf;
    subplot(2,1,1);
    imagesc(x,y,wn'); colorbar;caxis([-0.02 0.02])
    title(sprintf('t=%.2f',t));
    subplot(2,1,2);
    mesh(x,y,wn'); colorbar;caxis([-0.02 0.02])
    axis([0 Lx 0 Ly -0.05 0.05]);
    shg; pause(0.01);
end