clear all
close all
clc

dx = 0.1
dt=0.001
X1=0: dx: 1
X2=0: dx: 1
T = 0: dt: 1
M=zeros(length(X1),length(X2))
for i=1:1:length(X1)
    for j=1:1:length(X2)
        M(i,j)= sin(pi*X1(i))*cos(pi*X2(j))+1;
    end
end
lista = {M}
for k=2: 1: length(T)
    a = lista{end};
    b = zeros(size(a));
    for i=2:1:length(X1)-1
        for j=2:1:length(X2)-1
            b(i,j) = dt/(dx)^2 *(a(i-1,j)+a(i+1,j)+a(i,j-1)+a(i,j+1)-4*a(i,j))+a(i,j);
        end
    end
    b(:,1) = b(:,2);
    b(:,end) = b(:,end-1);
    b(1,:) = b(2,:);
    b(end,:) = b(end-1,:);
    lista{k} = b;
end
for t=1:1:length(lista)
    mesh(lista{t})
    zlim([-1,2])
    drawnow
end