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 < 1/2 )
  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)   % Wire-mesh plot of the solution
  mesh(ttplot,[-70 30]);
  xlabel('Position');  ylabel('Time');
  title('Diffusion of a delta spike');
subplot(122)   % Contour plot of the solution
  cs = contour(rot90(ttplot),0:.5:5,xplot,tplot);
  xlabel('Position');  ylabel('Time'); clabel(cs,0:5);
  title('Contour plot');
subplot(111)