function   [T,W]=rksys57(yprime,t0,tlast,w0,N);
%%function   [T,W]=rksys57(yprime,t0,tlast,w0,N);
%%   Runge-Kutta for systems of first order initial value problems
%%                 uj'=fj(t,u1,...,um)   j=1,...,m
%%                 on [a,b] with uj(a)=w0(j), j=1,...,m;
%%   It is assumed that the functions fj are given in an appropriate
%%   m-file that produces derivatives when called upon.
%%   INPUTS
%% yprime - String containing name of user-supplied problem description.
%%          yprime = 'fun', 'fun' corresponds to a function file,
%%                   fn_name.m where z=fn_name(t,U),
%%          t      - Time (scalar).
%%          U      - Solution row-vector [u(1),u(2),...,u(m)];
%%          yprime - Returns derivative values [u'(1),u'(2),...,u'(m)]
%%                   as row vector from t and [u(1),...,u(m)] values
%% t0     - Initial value of t.
%% tlast  - Last value of t.
%% w0     - Initial value row-vector.
%% N      - Number of steps.
%% 
%% OUTPUT
%%   T    -  column vector of times
%%   W    -  matrix of solutions  ... each row corresponding to the
%%           appropriate time, each column to a uj

  h=(tlast-t0)/N;
  T(1,1)=t0;
  W=w0;

   for j=1:N;
      k1=h*feval(yprime,T(j,1),W(j,:));
      k2=h*feval(yprime,T(j,1)+h/2,W(j,:)+(k1/2));
      k3=h*feval(yprime,T(j,1)+h/2,W(j,:)+(k2/2));
      k4=h*feval(yprime,T(j,1)+h,W(j,:)+k3);
      W(j+1,:)=W(j,:)+(k1+2*k2+2*k3+k4)/6;
      T(j+1,1)=T(j,1)+h;
   end;
