% This program solves the heat equation on a rectangle 
% PDE: ut=k(uxx+uyy)+f(x,y,t) 0<x<1, 0<y<1
% BC: u=0 on all 4 sides
% IC: u(x,y,0)=2sin(pi*x)sin(pi*y)
% du(x,y,t)/dt= d^2 u(x,y,t)/dx^2 + d^2 u(x,y,t)/dy^2 
% + f(x,y,t) on the square 0 < x,y < 1.  
% A finite-difference scheme is used with mesh width dx=dy.
% Modifed from Edward Allen, TTU
function zz=heat2d;
clear
clf
tf=1; %final time
n=40; % grid is 40X40
dx=1/n; 
dy=dx;
kk=.1; % Thermal diffusivity
nt=4*tf*n*n;  
dt=tf/nt; %For stability, we need dt/(dx*dx)+ dt/(dy*dy) <= 1/2.
x=linspace(0,1,n+1);
y=linspace(0,1,n+1);
% The initial values u(x,y,0) are specified next.
u=2*sin(pi*x)'*sin(pi*y); %[u(x,y,t)=2*exp(-(2 pi^2 k t))*sin(pi*x)*sin(pi*y)]
nm=1+(n/2); % center of rectangular plate
 %Boundary conditions are specified next. 
 for i=1:n+1
   for j=1:n+1
      u(1,j)=0;
      u(i,1)=0;
      u(n+1,j)=0;
      u(i,n+1)=0;
   end
end
%Finite Difference Method
v=u;
t=0;
jj=0;
for k=1:nt
   jj=jj+1;
   for i=2:n
      for j=2:n
         hlp=u(i+1,j)-2*u(i,j)+u(i-1,j); %u_xx approximation 
         hlp=hlp+u(i,j+1)-2*u(i,j)+u(i,j-1); %u_xx+u_yy approximation
         v(i,j)=u(i,j) + dt*kk*hlp/(dx*dx) + dt*f(x(i),y(j),t); 
      end
   end
   t=t+dt;
   u=v; 
   khlp=n*floor(k/n);
   if(khlp==k)  
      surf(y,x,u);
      xlabel('y')
      ylabel('x')
      zlabel('u(x,y,t)')
      title(sprintf('t=%.2f', t));
      colorbar
      caxis([-2,2])
      axis([0 1 0 1 -2 2])
      pause(.01);
   end  
 end
 
 disp(sprintf(' %9.3f  %9.3f  %9.5f ',x(nm),y(nm),u(nm,nm))); 
 % External force may be time-depenent
%------------------------------------
function ff=f(x,y,t)
%ff=exp(3*t)*sin(3*pi*x)'*sin(5*pi*y); 
%ff= 6*sin(2*pi*(x+y+t));
ff=0;
return
