function pdex

% A code to solve the diffusion equation with a sawtooth initial condition

L=pi;
N = 100;
Dx = L/(N+1);
for n = 2:N-1
  D2(n,n)  = -2; D2(n-1,n) = 1; D2(n+1,n) = 1;
end
D2(1,1)=-2; D2(2,1)=1; D2(N,N)=-2; D2(N-1,N)=1;
D2 = 1/Dx^2*D2;

tt = 0.2; Nt = 1001;  Dt = tt/Nt;  
x = linspace(0,L,N); t = linspace(0,tt,Nt);  % Space/time grid

% The discretized diffusion term can be written as a matrix
% product, D2*u, to compactify the coding

u(:,1)= (mod(x*2-pi,pi*2)-pi)/2;    % Initial condition
for n=1:Nt
  u(:,n+1) = u(:,n)+D2*u(:,n)*Dt;   % Time step the PDE
end

ns=[1 11 101 401 1001];             % Some sample snapshots

subplot('position',[0.1 0.5 0.8 0.3])    % Ouput plots hereon
plot(x,u(:,ns),'r');
ylabel('u(x,t)');
axis([0 pi -2 2])
text(2,1,'t=0, 0.002, 0.02, 0.08 and 0.2')
text(.1,-1,'Numerical solution')

x=linspace(0,L,1000);
subplot('position',[0.1 0.15 0.8 0.3])
for nt=ns
  u=0*x;
  for n=1:100
    u = u-2*cos(n*pi/2)/n*sin(n*x)*exp(-n^2*t(nt));   % Fourier series solution
  end
  plot(x,u),hold on
end
hold off
axis([0 pi -2 2])
xlabel('Position x');
ylabel('u(x,t)');
text(.1,-1,'Fourier series solution with 100 terms')

