function txy = classRK4(fcn,tstart,tend,xstart,nsteps)
% usage: txy = classRK4(fcn,tstart,tend,xstart,nsteps)
% description: This routine uses the classical RK4 on two
% dimensional problems with stepsize h=(tend-tstart)/nsteps.
% Inputs are the function name f for the function f(t,x) to be
% used (x is a column 2-vector), start time tstart, end time
% tend, initial value vector xstart and number of steps nsteps.
% Output is (nsteps+1)x3 matrix txy of time, x, y values as rows,
% where x, y are phase plane variables, equivalent to x(1),x(2).

% local variables:
% h: stepsize
% t: current time
% k1, k2, k3, k4: intermediate vectors

% initialization
txy = zeros(3,nsteps+1); % so columns of txy are [time;xcoord;ycoord]
txy(1,:)=linspace(tstart,tend,nsteps+1);
txy(2,1) = xstart(1);
txy(3,1) = xstart(2);
h = (tend-tstart)/nsteps;
% main loop
for ii=1:nsteps 
    k1 = h*feval(fcn,txy(1,ii),txy(2:3,ii));
    k2 = h*feval(fcn,txy(1,ii)+h/2,txy(2:3,ii)+0.5*k1);
    k3 = h*feval(fcn,txy(1,ii)+h/2,txy(2:3,ii)+0.5*k2);
    k4 = h*feval(fcn,txy(1,ii)+h,txy(2:3,ii)+k3);
    txy(2:3,ii+1)=txy(2:3,ii)+(k1+2*k2+2*k3+k4)/6;
end
% finally, transpose the result to conform to Matlab conventions
txy = txy';