dftcs.m

% dftcs - Program to solve the diffusion equation 
% using the FTCS scheme
clear; help dftcs;  % Clear memory and print header
tau = input('Enter time step - ');
N = input('Enter the number of grid points - ');
L = 1.;  % The system extends from x=-L/2 to x=L/2
h = L/(N-1);  % Grid size
kappa = 1.;  % Diffusion coefficient
coeff = kappa*tau/h^2;
if( coeff < .5 )
  disp('Solution is expected to be stable');
else
  disp('WARNING: Solution is expected to be unstable');
end
%%%%% Initial conditions and boundary conditions %%%%%
tt = zeros(N,1); % Initialize temperature to zero at all points
tt_new = zeros(N,1);  % Initialize temporary array used by FTCS
tt(N/2) = 1/h;   % Initial cond. is delta function in center
%% The boundary conditions are tt(1) = tt(N) = 0
%%%%% Set up loop and plot variables %%%%%
xplot = (0:N-1)*h - L/2;   % Record the x scale for plots
iplot = 1;                 % Counter used to count plots
nstep = 250;               % Maximum number of iterations
nplots = 25;               % Number of snapshots (plots) to take
plot_step = nstep/nplots;  % Number of time steps between plots
%%%%% MAIN LOOP %%%%%
for istep=1:nstep
  % FTCS scheme
  tt_new(2:(N-1)) = tt(2:(N-1)) + ...
      coeff*(tt(3:N) + tt(1:(N-2)) - 2*tt(2:(N-1)));
  tt = tt_new;        % Reset temperature to new values
  if( rem(istep,plot_step) < 1 )   % Every plot_step steps
    ttplot(:,iplot) = tt(:);       % record tt(i) for plotting
    tplot(iplot) = istep*tau;      % Record time for plots
    iplot = iplot+1;
    fprintf('%g out of %g steps completed\n',istep,nstep);
  end
end
subplot(121)
  mesh(tplot,xplot,ttplot);
  view([-70 30]);
  xlabel('Time');  ylabel('x');
  title('Diffusion of a delta spike');
subplot(122)
  cs = contour(xplot,tplot,flipud(rot90(ttplot)),0:.5:5);
  xlabel('x');  ylabel('Time'); clabel(cs,0:5);
  title('Contour plot');
subplot(111)