function [time, zzrk4, nfeval] = rk4(zzero,tfinal,nsteps)
% Runge-Kutta integrator for n-dimensional system of ODEs
% [time, zzrk4, nfeval] = rk4(zzero,tfinal,nsteps);
% input parameters
%   zzero    starting vector (column vector)
%   tfinal   end of time interval 
%   nsteps   number of time steps
% output parameters
%   time     vector of discrete times
%   zzrk4     matrix of solution values
%   nfeval   number of function evaluations
n = length(zzero);
fprintf('RK4 iteration in %g dimensions ...',n)
zz=zzero; t=0;
itntable(n,t,zz)
time=[t]; zzrk4=[zz]; nfeval=0;
h=tfinal/nsteps;
% time step loop
for k=1:nsteps;
K1 = rhsf(t,zz); nfeval=nfeval+1;
K2 = rhsf(t+0.5*h,zz+0.5*h*K1); nfeval=nfeval+1;
K3 = rhsf(t+0.5*h,zz+0.5*h*K2); nfeval=nfeval+1;
K4 = rhsf(t+h,zz+h*K3); nfeval=nfeval+1;
    t=t+h; zz = zz + (h/6)*(K1+2*K2+2*K3+K4);
    time=[time,t]; zzrk4=[zzrk4,zz];
    itntable(n,t,zz)
end
% end of time steps
fprintf('all done\n')
return

function itntable(n,t,zz)
if n==2,  % two-dimensional problem
   if t==0; % print header
   fprintf('\n  t        x           y          \n')
   fprintf(' %5.2f  %10.6f  %10.6f   \n',t,zz(1),zz(2))
   end
   fprintf(' %5.2f  %10.6f  %10.6f   \n',t,zz(1),zz(2))
end
return
