%% 2D BVP / Parabolic PDE examples - Converging to Steady-State (Heat Eqn.)
% Solving Heat equation dT/dt=alpha(Txx+Tyy) subject to boundary conditions (9.1)

%Domain in Centimeters
w=10;       %Width
h=15;       %Height
T0=100;     %Temperature Peak
dx=2.5;    
dy=dx;
x=[0:dx:w]; 
y=[0:dy:h];
[X,Y]=meshgrid(x,y);
Ta=T0.*sinh(pi*Y./w).*sin(pi*X./w)/sinh(pi*h/w)


%% Apply 2D PDE solution, in time

n=length(x);
m=length(y);
T1=zeros(size(X));
T1(end,:)=T0*sin(pi*x/(w));
T=zeros(size(X));
toler=1e-2;
t=0;
kappa=0.5;  %thermal diffusivity of copper in cm^2/s
dt=0.5*min([dx,dy])^2/(2*kappa)
while abs(norm(T1)-norm(T))>toler
    T=T1;
    t=t+dt;
    for i=2:m-1
        for j=2:n-1
            T1(i,j)=T(i,j)+(kappa*dt/(dx^2))*(T(i-1,j)+T(i+1,j)-2*T(i,j))+(kappa*dt/(dy^2))*(T(i,j-1)+T(i,j+1)-2*T(i,j));
        end;
    end;    
    %Plot Results
    figure(3);
    imagesc(x,y,T1); 
    axis xy;
    xlabel('Width (cm)');
    ylabel('Height (cm');
    title(['Forward Euler Numerical Solution: Temperature (C), t=',num2str(t)]);
    colorbar;
    pause(0.01);
end;

T1
t
t/dt