% Solution to 2D Heat Equation u_xx+u_yy=u_t, u(x,y,0)=0, Animated
% u(0,y,t)=u(pi,y,t)=u(x,0,t)=(x,pi,t)=0
% David Bloore, August 2011

spd=input('Speed (1,2,3,4)?');

switch spd
    case 1
        max=13; % highest m,n --fast 13 --slow 25
        d=80; % number of evaluation points along each axis --fast 80 --slow 300
        i=150; % number of intervals of time --fast 150 --slow 600
        tend=2.5; % end time (start time=0)
        dt=tend/i; % incriment of time
    case 2
        max=17;
        d=120;
        i=200;
        tend=1.5;
        dt=tend/i;
        
    case 3
        max=21;
        d=200;
        i=300;
        tend=1.5;
        dt=tend/i;
    case 4
        max=25;
        d=300;
        i=600;
        tend=1.5;
        dt=tend/i;
end

[X,Y]=meshgrid(linspace(0,pi,d),linspace(0,pi,d));
Z=zeros(d,d);

subplot('Position',[0.05 0.1 0.85 0.8]);
ss=surf(X,Y,Z);
title('Solution to 2D Heat Equation');
xlabel('0 \leq x \leq \pi');
ylabel('0 \leq y \leq \pi');
zlabel('u(x,y)');
set(ss,'EdgeColor','none');
xlim([0 pi]);
ylim([0 pi]);
zlim([0 1.4]);
colorbar('Location','Eastoutside');
set(gcf,'Position',[0 356 1440 900]);
caxis([0 1.4]);
light('Position',[0 -6 6]);

tpx=0;t=0;
subplot('Position',[0.94 0.1 0.02 0.8]);
title('Time');
tp=plot(tpx,t,'o');
xlim([-0.1 0.1]);
ylim([-0.1 tend+0.1]);
grid on;

for t=0:dt:tend;
    for mm=1:2:max;
        Z=zeros(d,d);
        for n=1:2:mm;
            for m=1:2:mm;
                Z=Z+(exp(-m^2*t)*exp(-n^2*t))*(16/pi^2)*((sin(m*X)/m).*(sin(n*Y)/n));
            end
        end
    end
    set(ss,'ZData',Z);
    set(tp,'YData',t);
    drawnow;
end