% 2-D wave equation

% 4x4 grid with Dirichlet BC and Neumann BC

nx=4;
ny=4;

un = zeros(nx,ny);     % array of displacement values
unm1 = zeros(nx,ny);
unp1 = zeros(nx,ny);

nt = 31;

g=zeros(nt+1,1);

dx=1;       % grid spacing
dy=1;

dt = .2;

T = 1; d = 1;

C = (T/d)*dt^2/dx^2;

f = [0 0 0 0; 0 10 10 0; 0 10 10 0; 0 0 0 0]; % forcing fcn



figure(1)
mesh(un)

for n=1:nt
   for i = 2:nx-1     % do explicit calculation
       for j = 2:nx-1
          unp1(i,j) = 2*un(i,j) - unm1(i,j) ...
          + C*(un(i+1,j) + un(i-1,j) ...
          + un(i,j+1) +  un(i,j-1)- 4*un(i,j)  ) ...
          - f(i,j)*dt^2/1;
       end
   end
   
   unp1(1,2) = unp1(2,2);
   unp1(1,3) = unp1(2,3);
    
    unm1=un; % replace unm1 with un
    un=unp1; % replace un with unp1

    figure(n+1)
    g(n+1)=un(2,3);
    mesh(un) 
  
end
figure(n+2)
plot(0:dx:nt,g)