%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Backward Euler method for the solution of ODEs of the form:
% y'(t)  =  f(t,y(t))
% The entries of this function are:
% a) Initial time
% b) Final time
% c) Number of time steps between initial and final time steps
% d) The function f(t,y(t))
% e) The gradient of f(t,y(t)) with respect to y(t)
% e) Initial condition y(t=0)=y0
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
function [y,t]  =  CrankNicholson(t0,tf,NTime,f,df,y0)

Dt  =  (tf-t0)/(NTime-1);  %  Time step constante
t   =  (linspace(t0,tf,NTime))';  %  Variable independiente, normalmente tiempo

NVariables  =  length(y0);
y           =  zeros(NTime,NVariables);             %  Inicialización de la variable dependiente


y(1,:)      =  y0;

g     =  @(y,yn,t,Dt)  y - yn - (Dt/2)*(f(t,y) + f(t-Dt,yn));
if NVariables==1
   DgDy  =  @(y,t,Dt)  1 - (Dt/2)*df(t,y);
else
   DgDy  =  @(y,t,Dt)  eye(length(y)) - (Dt/2)*df(t,y);    
end
tol   =  1e-6;
NMaxIter  =  100;
for itime=2:NTime
    y(itime,:)  =  MyNewtonRaphsonODE(g,DgDy,tol,NMaxIter,y(itime-1,:)',t(itime),Dt);
    
    
end

 
