function [t,y] = odeSchemes(odefun,tlim,y0,Nt,scheme)
% odeSchemes.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

% Time discretisation
t0 = tlim(1);
tf = tlim(2);
t  = (t0:(tf-t0)/Nt:tf)';

% Initialisation
y      = zeros(Nt+1,length(y0));
y(1,:) = y0;
h      = t(2)-t(1);

%%% Forward euler
if strcmp(scheme,'euler')
    for j = 2:Nt+1
        y(j,:) = y(j-1,:) + h * feval(odefun,t(j-1),y(j-1,:));
    end
  
%%% Midpoint
elseif strcmp(scheme,'midpoint')    
    for j = 2:Nt+1
        k1     = feval(odefun,t(j-1),y(j-1,:));
        k2     = feval(odefun,t(j-1)+h/2,y(j-1,:)+h/2*k1);
        y(j,:) = y(j-1,:) + h * k2;
    end
    
%%% Heun's
elseif strcmp(scheme,'heuns')
    for j = 2:Nt+1
        k1     = feval(odefun,t(j-1),y(j-1,:));
        y1     = y(j-1,:) + h * k1;
        k2     = feval(odefun,t(j-1)+h,y1);
        y(j,:) = y(j-1,:) + h * (k1+k2)/2;
    end
    
%%% RK4
elseif strcmp(scheme,'RK4')
    for j = 2:Nt+1
        k1     = feval(odefun,t(j-1),y(j-1,:));
        k2     = feval(odefun,t(j-1)+h/2,y(j-1,:)+h/2*k1);
        k3     = feval(odefun,t(j-1)+h/2,y(j-1,:)+h/2*k2);
        k4     = feval(odefun,t(j-1)+h  ,y(j-1,:)+h*k3);
        y(j,:) = y(j-1,:) + h * (1/6*k1 + 1/3*k2 + 1/3*k3 + 1/6*k4);
    end
    
%%% Unkown scheme
else
    error('odeSchemes.m : unknown scheme')
end
end
