% A. Donev, PDE Spring 2016
% The problem we are solving is the heat equation with Dirichlet boundary conditions
% The same as in the Maple worksheet:
% PDE: diff(u(x, t), t) = k*diff(u(x, t), x, x)
% BCs: u(0, t) = 0, u(pi, t) = 0
% IC: u(x,0) = piecewise(x<1/2,0,x>pi/2,pi,x=pi/2,pi/2);

% Recall the analytical solution is the infinite series:
% u(x,t) = sum((2*(cos((1/2)*pi*n)-(-1)^n))/n*sin(n*x)*exp(-n^2*t), n=1..infinity)
% In this Matlab code we will only use Matlab to evaluate and plot this solution

clear all;
N=100; % Choose how many terms to retain in the sum
M=1000; % How many points to use in the plot

u=zeros(M,1); % Allocate storage for the solution
lambdas=zeros(N,1); % Allocate memory for eigenvalues
coeffs=zeros(N,1); % and also for the Fourier coefficients

ns=[1:N]; % Vector of integers [1,2,...,N]

% Instead of something like
% for n=1:N; lambdas[n]=n^2; end;
% in Matlab it is much better (simpler, faster) to use dot-operator notation:
lambdas = ns.^2;
coeffs = (2*(cos((1/2)*pi.*ns)-(-1).^ns))./ns; % Fourier coefficients

% For the sum we will use a for loop
x=linspace(0,pi,M); % Vector of x values for the plotting grid

figure(1); clf
pause(5.0);

dt=0.001; % Time step for animation
for timestep=0:10
   t=timestep*dt;

   u=0;
   for n=1:N
     u = u + coeffs(n) * exp(-lambdas(n)*t) * sin(n*x);
   end

   plot(x, u, 'k-');
   xlim([0,pi]); ylim([-0.5,pi+0.5]);
   pause(0.5);
end
