clear
nx=250;
nt=350;
dt = .005;

dx = 2.5/(nx-1);
dy = 2.5/(nx-1);

u=zeros(nx,nx);
uo=zeros(nx,nx);
un=zeros(nx,nx);


c=1;
lc=1;
dlc=5;
R = 20;
for i=1:nx      % set up initial pulse
    for j=1:nx
        d = sqrt((i-125)^2 + (j-125)^2);
        if  d<R
            u(i,j) = .2*sin(pi/2+d/R*pi/2);
            uo(i,j) = u(i,j);
        end
    end
end
 
D=[0 1 0; 1 -4 1; 0 1 0]; % 2d Laplace operator

c1=dt^2*c^2/dx^2;

for n=1:nt
    
    if(1==1)     % use Laplace operator
        un=2*u - uo + c1*conv2(u,D,'same'); % new
    else
        
      for i = 2:nx-1     % do explicit calculation
          for j = 2:nx-1
             un(i,j) = 2*u(i,j) - uo(i,j) ...
             + c1*(u(i+1,j) -2*u(i,j) +  u(i-1,j) ...
              + u(i,j+1) -2*u(i,j) +  u(i,j-1) );
          end
       end
    end
        
    if mod(lc-1,dlc)==0     % animate graph
        figure(lc)
        mesh(u)   %plot(u(:,125))
        shg
    end
    lc=lc+1;
    
    uo=u; % curent become old
    u=un; % new become current
    
   
end