function [t,u]=cranknic(t0,T,N,f,y0)
%CRANKNIC Il metodo di Crank-Nicolson.
%   [T,U]=CRANKNIC(T0,T,N,F,Y0) risolve con 
%   il metodo di Crank-Nicolson il problema
%   di Cauchy y'(t)=F(t,y(t)) con y(T0)=Y0,
%   utilizzando N sotto-intervalli equispaziati.
%
%   Si vedano anche EULESP, EULIMP, PREDCOR.
%
%   Necessita del toolbox SYMBOLIC.

%   Saleri Fausto 24-01-01

h=(T-t0)/N; 
tt=[t0:h:T]; 
u(1)=y0;
syms x;  
t=t0; 
y = y0; 
fold = eval(f);
for t=tt(2:end)
  y = x;
  fun = ['x-',num2str(h,16),'*0.5*(',char(eval(f)),'+(',...
        num2str(fold,16),'))-(',num2str(u(end),16),')'];
  u = [u, zerof(fun,u(end))];
  y = u(end); fold = eval(f);
end

t=tt;

return
