N = 50;
map = zeros(N);

%map(5:8,5:8) = 1;

dt = 0.1;
dx = 1;

tempmap = map;
for t = 0:dt:50
    for i = 2:N-1
        for j = 2:N-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;
    mesh(map);
    zlim([0 1]);
    pause(0.05);
end
