% ode.m - Matthieu Aussal - Copyright 2018 - Licence CC BY-NC-SA 4.0
% CMAP - Ecole Polytechnique, Route de Saclay, 91128 Palaiseau CEDEX France
% For any informations, please contact  : matthieu.aussal@polytechnique.edu

clear all
close all
clc

%%%%%%%%%%%%%%%%%%%%%%%%% EXERCICE 1
% Define ODE
dty = @(t,y) -y;                        % Derivative
y0  = 1;                                % Initial condition

% Time interval
t0   = 0;
tf   = 2;

% Exact solution
Y = @(t) y0 * exp(-t);

% Analytic soution
Tref = linspace(t0,tf,1000);
Yref = Y(Tref);

% Numerical solution
N = 10;
tic
[Tnum,Ynum] = odeSchemes(dty,[t0 tf],y0,N,'euler');
toc

% Graphical representation
figure
plot(Tref,Yref)
hold on
plot(Tnum,Ynum,'o--')
xlabel('t')
ylabel('y(t)')
grid on
legend({'Analytic','Forward'})

% Errors
disp('Error max between forward Euler scheme and analytical solution :')
err = max(abs(Y(Tnum) - Ynum))


%%%%%%%%%%%%%%%%%%%%%%%%% Exercice 2a
% Operator
h = Tnum(2) - Tnum(1);
e = ones(N,1);
A = spdiags(e * [-1 1]/h,[-1 0],N,N);
B = spdiags(e,-1,N,N);

% Right hand side
F = zeros(N,1);

% Resolution
Ystar = (A+B)\F


%%%%%%%%%%%%%%%%%%%%%%%%% Exercice 2b
% Right hand side
F(1) = (1-h)/h*y0;

% Resolution
tic
Ystar = (A+B)\F;
toc

% Add initial condition
Yfwd  = [y0;Ystar];    

% Error
disp('Error max between matrix resolution and forward Euler scheme :')
err = max(abs(Yfwd - Ynum))


%%%%%%%%%%%%%%%%%%%%%%%%% Exercice 3
% Backward osperator
B = spdiags(e,0,N,N);

% Right hand side
F(1) = y0/h;

% Resolution
Ystar = (A+B)\F;

% Add initial condition
Ybwd  = [y0;Ystar];  

% Graphical representation
plot(Tnum,Ybwd,'+--')
grid on
legend({'Analytic','Forward','Backward'})

% Error
disp('Error max between matrix resolution and analytical solution :')
err = max(abs(Y(Tnum) - Ybwd))


%%%%%%%%%%%%%%%%%%%%%%%%% Exercice 4
% Centered scheme
A = spdiags(e * [-1 1]/(2*h),[-1 1],N,N);
B = spdiags(e,0,N,N);

% Right hand side
F(1) = y0/(2*h);

% Resolution
Ystar = (A+B)\F;

% Add initial condition
Yctr  = [y0;Ystar]; 

% Graphical representation
plot(Tnum,Yctr,'+--')
grid on
legend({'Analytic','Forward','Backward','Centered'})

% Error
disp('Error max between matrix resolution and analytical solution :')
err = max(abs(Y(Tnum) - Yctr))



