clear all
close all
%------------------------------
% Solution of the 2d wave equation by use of explicit 5-point method
%---------------------------------

L=1; % length of the string
c=1;
M=200; % amount of points on the (0,L)
x=linspace(0, L, M+1);
y=linspace(0, L, M+1);
dx=x(2)-x(1);% space grid size
dy=y(2)-y(1);
a=0.25;
dt=(a*dx)/c; % time grid size
T=400; % amount of time steps
tt=10;
s=2; % shift
rad=0.01;
ys=0.15;%23%1.1;
ob=2;
%-----------------------------
u=zeros(M+1,M+1,T+1);
[x,y]=meshgrid(0:dx:L);
f=inline('4*z.^2.*(1-z).*k.*(1-k)','z','k'); % an initial wave u(x,)
%f=inline('sin(2*pi*z).*cos(2*pi*z)','z','k'); % an initial wave u(x,)
u(:,:,1)=f(x,y);
u(1,1,1)=0; u(M+1,1,1)=0; u(M+1,M+1,1)=0; u(1,M+1,1)=0;

surf(x,y,f(x,y));
xlabel('x')
       ylabel('y')
       zlabel('u(x,y,t)')
shading interp
axis([0-rad L+rad 0-rad L+rad -ys ys])
% %-----------------------------
 %Calculate the second time row
 u(1,1,2)=0; u(M+1,1,2)=0; u(M+1,M+1,2)=0; u(1,M+1,2)=0;
 for i=2:M
     for j=2:M
     u(i,j,2)=(1-a^2-a^2)*u(i,j,1)+(a^2*(u(i+1,j,1)+u(i-1,j,1)))/2+(a^2*(u(i,j+1,1)+u(i,j+1,1)))/2;
     end;
 end;
 %-----------------------------
 for t=2:T
    u(1,1,t+1)=0; u(M+1,1,t+1)=0; u(M+1,M+1,t+1)=0; u(1,M+1,t+1)=0;
    for i=2:M
        for j=2:M
        u(i,j,t+1)=-u(i,j,t-1)+2*(1-a^2-a^2)*u(i,j,t)+a^2*(u(i+1,j,t)+u(i-1,j,t))+a^2*(u(i,j+1,t)+u(i,j-1,t));
        end;
    end;
    if mod(t,tt)==0
      surf(x,y,u(:,:,t+1));
      %view(8, 8);
%      colorbar
      shading interp  
      axis([0-rad L+rad 0-rad L+rad -ys ys])
      title(['t=' num2str(t)],'Color','b')
       xlabel('x')
       ylabel('y')
       zlabel('u(x,y,t)')
       hold on 
       pause(0.5)
       hold off
    end;
 end;
 %-----------------------------
 