Sample m-file

% diffuse1.m - Program to solve the diffusion equation
%              using the FTCS method.

clear;             % clear memory
help diffuse1.m    % print help information

%%%%%%%% INITIALISATION %%%%%%%%%%

% set up input variables
tmax = 1;                          % window size in t
nt = 100;                          % number of grid points in t
delta_t = tmax/nt;                 % step-size in t
nplots = 21;                       % number of plot points in t
nplotstep = nt/nplots;             % plot step-size in t

xmax = 2;                          % window size in x
nx = 11;                           % number of grid points in x
delta_x = xmax/(nx + 1);           % step-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 and zeros a column vector
A(floor((nx+1)/2)) = 1/delta_x;    % sets up an initial delta 'spike' in A
iplot = 1;                         % initialise 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)
    
    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 stuff
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 spike')% 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 spike')  % give the thing a title

subplot(111)                       % return plotting mode to one 
                                   % figure per window