%plot the analytic result with different number of modes
%problem4
close all

L=1;                %domain length
npts=100;           %discretization
x=0:L/(npts-1):L;   %space vector
t=[0,0.1,0.2,0.5,1];%time vector
N=[1,5,10,100];     %number of modes

D=0.5;  %diffusion constant
C0=1;   %initial concentration

A= @(n) 8*C0/((2*n-1)*pi)^3;
mode= @(n,t) sin((2*n-1)*pi/L*x)*exp(-D*t*((2*n-1)*pi/L)^2);

% figure(1)
% for i=1:length(N)       %different number of modes
%     subplot(2,2,i)
%     for j=1:length(t)   %loop over times
%         C=0;
%         for k=1:N(i)    %loop over modes
%             C=C + A(k)*mode(k,t(j));
%         end
%         plot(x,C)
%         hold on
%     end
% end

%----------------------
        
%compute 1 mode solution
N=1;
C1=zeros(length(t),npts);
for j=1:length(t)   %loop over times
    for k=1:N    %loop over modes
        C1(j,:)=C1(j,:) + A(k)*mode(k,t(j));
    end
end

%compute 2 mode solution
N=2;
C2=zeros(length(t),npts);
for j=1:length(t)   %loop over times
    for k=1:N    %loop over modes
        C2(j,:)=C2(j,:) + A(k)*mode(k,t(j));
    end
end

%compute 3 mode solution
N=3;
C3=zeros(length(t),npts);
for j=1:length(t)   %loop over times
    for k=1:N    %loop over modes
        C3(j,:)=C3(j,:) + A(k)*mode(k,t(j));
    end
end

figure(1)
subplot(2,3,1)
plot(x,C1)
title('1 mode solution (M1)')
xlabel('x')
ylabel('C(x,t)')
legend('0s','0.1s','0.2s','0.5s','1s')
set(gca,'YLim',[0.0,0.27])

subplot(2,3,2)
plot(x,C2)
title('2 mode solution (M2)')
xlabel('x')
ylabel('C(x,t)')
legend('0s','0.1s','0.2s','0.5s','1s')
set(gca,'YLim',[0.0,0.27])

subplot(2,3,3)
plot(x,C3)
title('3 mode solution (M3)')
xlabel('x')
ylabel('C(x,t)')
legend('0s','0.1s','0.2s','0.5s','1s')
set(gca,'YLim',[0.0,0.27])

subplot(2,3,4)
plot(x,C1-C2)
title('Difference M1 - M2')
xlabel('x')
ylabel('\Delta C(x,t)')
set(gca,'YLim',[-0.012,0.012])

subplot(2,3,5)
plot(x,C2-C3)
title('Difference M2 - M3')
xlabel('x')
ylabel('\Delta C(x,t)')

subplot(2,3,6)
plot(x,C1-C3)
title('Difference M1 - M3')
xlabel('x')
ylabel('\Delta C(x,t)')
set(gca,'YLim',[-0.012,0.012])
