% diffusion_01.m - Program to solve the diffusion equation
%                  using the Forward Time Centred Space (FTCS) method
%  The diffusion equation (LaTeX notation)  
%
%                 \partial_{t} f= \kappa \partial_{x}^2 f  
%
%  is solved by approximating it as:                    
%
%     f(x,t+dt)=f(x,t)+(\kappa*dt)( f(x+dx,t)-2 f(x,t)+f(x-dx,t))/dx^2
%
% Further information on the code is retrievable from
% http://www.physics.uq.edu.au/people/cochrane/phys3071/phys3071Notes/node8.html
    
  clc;               % clear screen
  clear all          % clear memory
    
%%%%%%%% INITIALIZATION %%%%%%%%%%
    
  tmax = 1;                          % horizon in t
  nt = 100;                          % number of grid points in t
  delta_t = tmax/nt;                 % mesh-size in t
  nplots = 21;                       % number of plots
  nplotstep = nt/nplots;             % number of steps before plotting 
  
  nx = 11;                           % number of grid points in x
  xmax = 2;                          % horizon in x
   
  delta_x = xmax/(nx + 1);           % mesh-size in x
   
  kappa = 1;                         % coupling constant
  c = kappa*delta_t/delta_x^2;       % constant of propagation
   
 % initial conditions and boundary conditions
  A = zeros(nx,1);                   % allocates zeros in a column vector
 
 % floor rounds a real to its integer part (by defect).
  A(floor((nx+1)/2)) = 1/delta_x;    % sets up an initial delta 'spike' in A
  iplot = 1;                       % initialize the plotting index
   
 % set up plotting matrices
 xplot = (1:nx)*delta_x - xmax/2;   % row vector of x-coordinates
 tplot(iplot) = 0;                  % assign first time-coordinate
 Aplot(:,iplot) = A(:);             % assign first value of A
    
    %%%%%%%% MAIN LOOP %%%%%%%%%%
    
% main loop
 for i = 1:nt
      
% evolve our solution using the FTCS method
    A(2:(end-1)) = A(2:(end-1)) + c*( A(1:(end-2)) + A(3:end) - 2*A(2:(end-1)) );
     
  % save current data for plotting
    if (rem(i,nplotstep) < 1)          % rem(a,b)=a-fix(a./b)*b
  %                                    % fix(x) rounds x to the integer
  %                                    % closest to zero
                                      
        iplot = iplot + 1;             % increment plotting index
        Aplot(:,iplot) = A(:);         % record current value of A for plotting
        tplot(iplot) = i*delta_t;      % record current time value
       
  % print some progress information
       fprintf('%g out of %g steps completed\n',i,nt);
      
    end
     
 end


  %%%%%%%% PLOTTING/OUTPUT %%%%%%%%%%
   
% plot 
 figure(1)                          % opens a new figure window
 subplot(121)                       % divides window into a 1x2 matrix 
                                      % and assigns next plot to the 
                                      % left-hand frame (frame 1)
 mesh(tplot,xplot,Aplot)            % wire mesh plot of the solution
 view([-70 30])                     % change the viewing direction
 xlabel('t')                        % label the x-axis
 ylabel('x')                        % label the y-axis
 title('Diffusion of a delta function: 2d plot')% give the thing a title
   
 subplot(122)                       % assigns next plot to right-hand 
                                      % frame (frame 2)
   
% find contours of 3D plot and assign to the vector cs
 cs = contour(xplot,tplot,flipud(rot90(Aplot)),0:0.1:2.5);
 xlabel('x')                        % label the x-axis
 ylabel('t')                        % label the y-axis
 clabel(cs,0:5)                     % label the contours
 title('Contour of a delta function')  % give the thing a title
   
   subplot(111)                       % return plotting mode to one 
                                      % figure per window
                                      
figure(2)                          % opens a new figure window
plot(xplot,Aplot(:,1:iplot))
xlabel('x')                           % label the x-axis
ylabel('f(x)')                        % label the y-axis
title('Diffusion of a delta function: plots at different times')% give the thing a title