%problem 5

close all
clear all

L=1;                %domain length
times=[0 0.1 0.2 0.5 1];%time vector
nt=length(times);
D=0.5;  %diffusion constant
C0=1;   %initial concentration

bc=zeros(nt,1);

% 3pts
npts=3;           %discretization
dx=L/(npts-1);
x3=0:dx:L;   %space vector
rhs = @(t,C) D*([0;C(1:end-1)] - 2.0*C + [C(2:end);0])./dx^2;
initC = C0*(x3(2:end-1)/L - x3(2:end-1).^2/L^2);

[t3,C3] = ode45(rhs,times,initC);
C3=[bc,C3,bc];

% 5pts
npts=5;           %discretization
dx=L/(npts-1);
x5=0:dx:L;   %space vector
rhs = @(t,C) D*([0;C(1:end-1)] - 2.0*C + [C(2:end);0])./dx^2;
initC = C0*(x5(2:end-1)/L - x5(2:end-1).^2/L^2);

[t5,C5] = ode45(rhs,times,initC);
C5=[bc,C5,bc];

% 11pts
npts=11;           %discretization
dx=L/(npts-1);
x11=0:dx:L;   %space vector
rhs = @(t,C) D*([0;C(1:end-1)] - 2.0*C + [C(2:end);0])./dx^2;
initC = C0*(x11(2:end-1)/L - x11(2:end-1).^2/L^2);

[t11,C11] = ode45(rhs,times,initC);
C11=[bc,C11,bc];

newPts=3*5*11;
xnew=0:1.0/(newPts-1):L;
C3new=zeros(length(times),newPts);
C5new=zeros(length(times),newPts);
C11new=zeros(length(times),newPts);
%difference
for i=1:length(times)
    C3new(i,:)=interp1(x3,C3(i,:),xnew);
    C5new(i,:)=interp1(x5,C5(i,:),xnew);
    C11new(i,:)=interp1(x11,C11(i,:),xnew);
end

%Plot
%-----
subplot(2,3,1)
plot(x3,C3)
title('3 points solution (P3)')
xlabel('x')
ylabel('C(x,t)')
legend('0s','0.1s','0.2s','0.5s','1s')

subplot(2,3,2)
plot(x5,C5')
title('5 points solution (P5)')
xlabel('x')
ylabel('C(x,t)')
legend('0s','0.1s','0.2s','0.5s','1s')

subplot(2,3,3)
plot(x11,C11)
title('11 points solution (P11)')
xlabel('x')
ylabel('C(x,t)')
legend('0s','0.1s','0.2s','0.5s','1s')


subplot(2,3,4)
plot(xnew,C3new-C5new)
title('Difference P3 - P5')
xlabel('x')
ylabel('\Delta C(x,t)')
set(gca,'YLim',[-0.07,0.02])

subplot(2,3,5)
plot(xnew,C5new-C11new)
title('Difference P5 - P11')
xlabel('x')
ylabel('\Delta C(x,t)')
set(gca,'YLim',[-0.07,0.02])

subplot(2,3,6)
plot(xnew,C3new-C11new)
title('Difference P3 - P11')
xlabel('x')
ylabel('\Delta C(x,t)')
set(gca,'YLim',[-0.07,0.02])
