clear all;
%% math 454 : plotting template
%
%   this script provides a template for plotting solutions of PDEs in one
%   space dimension using Fourier series.
%   The template solves the heat equation, but may be easily modified to
%   other PDEs whose solutions may be represented as Fourier series (esp. the wave equation).
%
%   author : Peter Bosler, University of Michigan Department of Mathematics
%   
%   Step 1:  change lines 18, 20, 21, 23 to match the given problem
%   Step 2:  change lines 32 and 33 to the correct the Fourier coefficients
%   Step 3:  change the eigenfunctions in line 50
%   Step 4:  push 'Run', save and print the figures
%
%% Change this section to match your problem

sumLimit = 40;  % number of terms in finite Fourier series 

L = 3;         % right endpoint of spatial domain 
k = 1;       % material constant 

t = [0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0];  % times to plot solutions

%% Fourier coefficients (change to match your problem)

a0 = 0;              % define the a0 coefficient here

a = zeros(sumLimit,1); % vector to hold cosine coefficients (don't change this line)
b = zeros(sumLimit,1); % vector to hold sine coefficiients  (don't change this line)
for n = 1 : sumLimit
   a(n) = 0;               % define the a(n) coefficient here
   b(n) = (-1)^(n-1)/n^2;      % define the b(n) coefficient here
end

%% spatial domain (not necessary to change, usually)

dx = 0.01;     % spacing of mesh points for discrete spatial domain
x = (0:dx:L)'; % spatial domain

%% calculate solutions (change eigenfunctions and time dependent coefficients)

u = zeros(length(x),length(t)) + a0;    % start with a0

for i = 1:length(t) % loop over time steps to plot
    
    for n = 1 : sumLimit % add the next term in the Fourier series
        
        % add appropriate eigenvalues & eigenfunctions to the next line
        u(:,i) = u(:,i) + a(n)*exp(-k*(n*pi/L)^2*t(i))*cos(n*pi*x/L) + b(n)*exp(-k*(n*pi/L)^2*t(i))*sin(n*pi*x/L);  
        
    end
end

%% plots

figure(1);clf;  % plot first 4 timesteps in figure 1
plot(x,u(:,1), x,u(:,2), x,u(:,3), x,u(:,4))
set(gca,'FontSize',16);
xlabel('x');
ylabel('u(x,t)');
ls1 = sprintf('t = %4.2f',t(1));
ls2 = sprintf('t = %4.2f',t(2));
ls3 = sprintf('t = %4.2f',t(3));
ls4 = sprintf('t = %4.2f',t(4));
legend(ls1,ls2,ls3,ls4);
xlim([0,L]);

figure(2);clf; % plot next 4 timesteps in figure 2
plot(x,u(:,5), x,u(:,6), x,u(:,7), x,u(:,8));
set(gca,'FontSize',16);
xlabel('x');
ylabel('u(x,t)');
ls1 = sprintf('t = %4.2f',t(5));
ls2 = sprintf('t = %4.2f',t(6));
ls3 = sprintf('t = %4.2f',t(7));
ls4 = sprintf('t = %4.2f',t(8));
legend(ls1,ls2,ls3,ls4);
xlim([0,L]);

figure(3);clf; % plot last 4 timesteps in figure 3
plot(x,u(:,9), x,u(:,10), x,u(:,11), x,u(:,12));
set(gca,'FontSize',16);
xlabel('x');
ylabel('u(x,t)');
ls1 = sprintf('t = %4.2f',t(9));
ls2 = sprintf('t = %4.2f',t(10));
ls3 = sprintf('t = %4.2f',t(11));
ls4 = sprintf('t = %4.2f',t(12));
legend(ls1,ls2,ls3,ls4);
xlim([0,L]);

