clear all;
close all;

N = 63*2+62*3;
M = 62*2+61*4;
map = zeros(N+2,M+2);
%map(2:N+1,2:M+1) = ones(N,M);
%map(2:(N+1)>N/2,2:(M+1)>M/2) = 1;
for i = 2:N+1
    for j = 2:M+1
        if (j>(1/50*((i-150).^2)+30)) && (j<(1/90*((i-150).^2)+190))
            map(i,j) = 1;
        end
    end
end

dx = 1;
dt = 0.1;

T = 0:dt:(100-dt);
tempmap = map;
tic
for t = 1:length(T)
    for i = 2:N+1
        for j = 2:M+1
            tempmap(i,j) = map(i,j)+dt*((map(i-1,j)-2*map(i,j)+map(i+1,j))/dx^2+(map(i,j-1)-2*map(i,j)+map(i,j+1))/dx^2);
        end
    end
    map = tempmap;
    %{
    map(rand(N)>0.999) = 1;
    map(5:10,5:8) = 1;
    map(5:6,9:12) = 1;
    map(9:10,9:12) = 1;
    mesh(map,ones(N+2,M+2));
    zlim([0 1]);
    imagesc(map)
    pause(0.1);
    %waitforbuttonpress;
    %}

end
toc
origmap = map(2:N+1,2:M+1);
mesh(origmap);
%%