function [t,u]=newmark(odefun,tspan,y0,Nh,param,...
                        varargin)
%NEWMARK  Risolve equazioni differenziali del II ord.
%  con il metodo di Newmark.
%  [T,Y]=NEWMARK(ODEFUN,TSPAN,Y0,NH,PARAM) con
%  TSPAN=[T0 TF] integra il sistema di equazioni dif-
%  ferenziali y''= f(t,y,y') dal tempo T0 a TF con
%  condizioni iniziali Y0=(y(t0),y'(t0)) utilizzando
%  il metodo di Newmark su una griglia equispaziata di
%  NH intervalli. Il vettore PARAM contiene, in ordine,
%  i parametri zeta e theta del metodo di Newmark.
%  La funzione ODEFUN(T,Y) deve ritornare uno scalare,
%  la variabile Y e' un array che contiene la funzione
%  soluzione in prima componente e la sua derivata in
%  seconda componente.
%  Ogni riga del vettore soluzione Y corrisponde ad
%  un istante temporale del vettore colonna T.
%  [T,Y] = NEWMARK(ODEFUN,TSPAN,Y0,NH,P1,P2,...) passa
%  i parametri addizionali P1,P2,... alla funzione
%  ODEFUN come ODEFUN(T,Y,P1,P2...).
tt=linspace(tspan(1),tspan(2),Nh+1);
y=y0(:); u=y.';
global glob_h glob_t glob_y glob_odefun;
global glob_zeta glob_theta glob_varargin glob_fn;
glob_h=(tspan(2)-tspan(1))/Nh;
glob_y=y; glob_odefun=odefun;
glob_zeta = param(1); glob_theta = param(2);
glob_varargin=varargin;
if ( ~exist( 'OCTAVE_VERSION' ) )
 options=optimset;
 options.Display='off';
 options.TolFun=1.e-12;
 options.MaxFunEvals=10000;
end
glob_fn =feval(odefun,tt(1),glob_y,varargin{:});
for glob_t=tt(2:end)
if ( exist( 'OCTAVE_VERSION' ) )
  w = fsolve('newmarkfun', glob_y );
else
  w = fsolve(@(w) newmarkfun(w),glob_y,options);
end
  glob_fn =feval(odefun,glob_t,w,varargin{:});
  u = [u; w.']; glob_y = w;
end
t=tt;     clear glob_h glob_t glob_y glob_odefun;
clear glob_zeta glob_theta glob_varargin glob_fn;
end

function z=newmarkfun(w)
 global glob_h glob_t glob_y glob_odefun;
 global glob_zeta glob_theta glob_varargin glob_fn;
 fn1=feval(glob_odefun,glob_t,w,glob_varargin{:});
 z(1)=w(1) - glob_y(1) -glob_h*glob_y(2)-...
   glob_h^2*(glob_zeta*fn1+(0.5-glob_zeta)*glob_fn);
 z(2)=w(2) - glob_y(2) -...
   glob_h*((1-glob_theta)*glob_fn+glob_theta*fn1);
end