% Solution of the 2D Heat Equation Using a Forward Difference Scheme
% u_t = k (u_xx + u_yy)

% Initialize Data
%      L, W   = Length and Width of Playe, 
%      T      = for Time Interval [0, T]
%      Nx, Ny = Number of Points in Space Grid, 
%      M      = Number of Time Steps
%      dx, dy = Delta x and Delta y.
% Set your own values of L, W, T, k 
clear
L=1;
W=2;
T=1;
k=1/10;
Nx=20;
Ny=20;
dx=L/Nx;
dy=W/Ny;

% Find stable timestep
dt=min([dx,dy])^2/k/4;
M=ceil(T/dt);
alphax = k*dt/dx^2;
alphay = k*dt/dy^2;

% Spatial grid
[x,y]=meshgrid(0:dx:L,0:dy:W);

% Initialize u0 and u1, Change to your initial condition
u0=sin(pi*x/L).*sin(pi*y/W);         % Fixed conditions all around
% u0=y.*(2-y).*(1-x.^2);             % Example using Neumann at x=0  
u1=zeros(Nx+1,Ny+1);
H=max(max(u0));

% Plot initial condition
surf(x,y,u0,'FaceColor','red','EdgeColor','none')
camlight left;
lighting phong
xlabel('x')
ylabel('y')
title('Solution at t = 0')
axis([0,L,0,W,0,H])
frame=1;
Mov(frame)=getframe(gcf); 
pause (1)
 
% Partial Difference Equation (Numerical Scheme)

for m=1:M
   for i=2:Nx 
       for j=2:Ny
           u1(i,j)=u0(i,j)+alphax*(u0(i+1,j)-2*u0(i,j)+u0(i-1,j)) ...
               +alphay*(u0(i,j+1)-2*u0(i,j)+u0(i,j-1));
       end
   end
   % Boundary conditions - Set your own
   % Example: For thermally insulated side at x=0, replace  u1(1,:)=0;
   % with  u1(1,:)= u1(2,:);
   u1(1,:)=0;           % Condition at y=0
   u1(Nx+1,:)=0;        % Condition at y=W
   u1(1,:)=0;           % Fixed Condition at x=0 for 1st example IC
   %u1(:,1)=u1(:,2);    % Insulating Condition at x=0 for 2nd example IC
   u1(:,Ny+1)=0;        % Condition at x=L
   u0=u1; 
   
   % Plot 3D solution every 10 time steps
   if mod(m,10)==0
       frame=frame+1;
       surf(x,y,u1,'FaceColor','red','EdgeColor','none')
       camlight left;
       lighting phong
       xlabel('x')
       ylabel('y')
       title(['Solution at t = ' num2str(m*dt)])
       axis([0,L,0,W,0,H])
       Mov(frame)=getframe(gcf); 
       pause (0.5)
   end
end

% Extra - show or create a movie
%    movie(gcf, Mov)         % plays movie
%    movie(gcf, Mov,1,2)     % plays at 2 fps 
%    movie(gcf, Mov,10,5)    % repeats 10 times at 5 fps
% Create movie 
%    movie2avi(Mov, 'heat2d.avi', 'compression', 'None');