function [t,u]=cranknic(odefun,tspan,y,Nh,varargin)
%CRANKNIC  Solve differential equations using the
%  Crank-Nicolson method.
%  [T,Y]=CRANKNIC(ODEFUN,TSPAN,Y0,NH) with TSPAN=[T0,TF]
%  integrates the system of differential equations
%  y'=f(t,y) from time T0 to TF with initial condition
%  Y0 using the Crank-Nicolson method on an equispaced
%  grid of NH intervals.Function ODEFUN(T,Y) must return
%  a column vector corresponding to f(t,y). Each row in
%  the solution array Y corresponds to a time returned
%  in the column vector T.
%  [T,Y] = CRANKNIC(ODEFUN,TSPAN,Y0,NH,P1,P2,...) passes
%  the additional parameters P1,P2,... to the function
%  ODEFUN as ODEFUN(T,Y,P1,P2...).
tt=linspace(tspan(1),tspan(2),Nh+1);
u(:,1)=y;
global the_h the_t the_y the_odefun;
the_h=(tspan(2)-tspan(1))/Nh;
the_y=y;
the_odefun=odefun;

fid=fopen('myfun.m','w');
fprintf(fid,'function z=myfun(w)\n');
fprintf(fid,'global the_h the_t the_y the_odefun;\n');
fprintf(fid,['z=w-the_y-0.5*the_h*(feval(the_odefun,the_t,w)+',...
    'feval(the_odefun,the_t,the_y));\n']);
fprintf(fid,'return\n');
fclose(fid);
if( ~exist('OCTAVE_VERSION') )
    options=optimset('fsolve');
    options.Display='off';
    options.TolFun=1.e-06;
    options.MaxFunEvals=10000;
end
for the_t=tt(2:end)
    if ( exist('OCTAVE_VERSION') )
      [w, info] = fsolve("myfun",the_y);
    else
        w = fsolve('myfun',the_y,options);
    end
    u = [u w];
    the_y = w;
end
t=tt;
clear the_h the_t the_y the_odefun;
return