function [t,u]=cranknic(odefun,tspan,y,Nh,varargin)
%CRANKNIC    Risolve equazioni differenziali usando il metodo di Crank-Nicolson.
%   [T,Y] = CRANKNIC(ODEFUN,TSPAN,Y0,NH) with TSPAN = [T0 TFINAL] 
%   integra il sistema di equazioni differenziali y' = f(t,y) dal  
%   tempo T0 a TFINAL con condizione iniziale Y0 usando il metodo di Crank-
%   Nicolson su una griglia di NH intervalli equispaziati.
%   La funzione ODEFUN(T,Y) deve ritornare un vettore colonna corrispondente
%   a f(t,y). Ogni riga dell'array soluzione Y corrispe alla soluzione calcolata
%   per il corrispondente vettore colonna T.
%   [T,Y] = CRANKNIC(ODEFUN,TSPAN,Y0,NH,P1,P2,...) passa i parametri addizionali
%   P1,P2,... alla funzione ODEFUN come ODEFUN(T,Y,P1,P2...). 

%   A.Quarteroni,F.Saleri, Introduzione al Calcolo Scientifico, 2004

h=(tspan(2)-tspan(1))/Nh;
tt=linspace(tspan(1),tspan(2),Nh+1);
u(1)=y; 
fold = feval(odefun,tspan(1),y,varargin{:});
syms x;  
y = x;
for t=tt(2:end)
    y = x;
    fun = inline(['x-',num2str(h,16),'*0.5*(',char(feval(odefun,t,y,varargin{:})),'+(',...
            num2str(fold,16),'))-(',num2str(u(end),16),')'],'x');
    u = [u, fzero(fun,u(end))];
    fold = feval(odefun,t,u(end),varargin{:});
end
t=tt;

return
