L = 1;
N = 30;
h = L/N;
t_f = .06;
M = 250;
k = t_f/M;

x = 0:h:L;
y = x;
% define the initial data:
K1 = 2;
K2 = 4;
for i=0:N
    for j=0:N
        u0(i+1,j+1) = sin(K1*pi*x(i+1)/L)*sin(K2*pi*y(j+1)/L);
    end
end
[u,t] = two_d_heat(L,t_f,h,k,u0);
figure(1)
plot(t,squeeze(max(max(u))))
figure(2)
plot(t,log10(squeeze(max(max(u)))))
figure(3)
surf(u(:,:,1))
figure(4)
surf(u(:,:,101));
figure(5)
surf(u(:,:,251));

% repeat with coarser time-step
M = 20;
k = t_f/M;
[u,t] = two_d_heat(L,t_f,h,k,u0);
figure(1)
plot(t,squeeze(max(max(u))))
figure(2)
plot(t,log10(squeeze(max(max(u)))))
figure(3)
surf(u(:,:,1))
figure(4)
surf(u(:,:,11));
figure(5)
surf(u(:,:,21));