function [tvals,yvals]=FixedRK(fname,t0,y0,h,k,n)
tvals=[];yvals=[];
// 
// Produces approximate solution to the initial value problem
// 
//      y'(t) = f(t,y(t))     y(t0) = y0
// 
// using a strategy that is based upon a k-th order
// Runge-Kutta method. Stepsize is fixed.
// 
// Pre:  fname = function f.
//       t0 = initial time.
//       y0 = initial condition vector.
//       h  = stepsize.
//       k  = order of method. (1<=k<=5).
//       n  = number of steps to be taken,
// 
// Post: tvals(j) = t0 + (j-1)h, j=1:n+1
//       yvals(:j) = approximate solution at t = tvals(j), j=1:n+1
// 
tc = t0;
yc = y0;
tvals = tc;
yvals = yc;
fc = fname(tc,yc);
for j = 1:n
  [tc,yc,fc] = RKstep(fname,tc,yc,fc,h,k);
  yvals = [yvals,yc];
  tvals = [tvals,tc];
end


function [tnew,ynew,fnew]=RKstep(fname,tc,yc,fc,h,k)
tnew=[];ynew=[];fnew=[];
// 
// Pre:  fname is a function of the form f(t,y)
//       where t is a scalar and y is a column d-vector.
// 
//       yc is an approximate solution to y'(t) = f(t,y(t)) at t=tc.
// 
//       fc = f(tc,yc).
// 
//       h is the time step.
// 
//       k is the order of the Runge-Kutta method used, 1<=k<=5.
// 
// Post: tnew=tc+h, ynew is an approximate solution at t=tnew, and
//       fnew = f(tnew,ynew).
 
if k==1 then
  k1 = h*fc;
  ynew = yc+k1;
   
elseif k==2 then
  k1 = h*fc;
  k2 = h*fname(tc+h,yc+k1);
  ynew = yc+(k1+k2)/2;
   
elseif k==3 then
  k1 = h*fc;
  k2 = h*fname(tc+h/2,yc+k1/2);
  k3 = h*fname(tc+h,yc-k1+2*k2);
  ynew = yc+(k1+4*k2+k3)/6;
   
elseif k==4 then
  k1 = h*fc;
  k2 = h*fname(tc+h/2,yc+k1/2);
  k3 = h*fname(tc+h/2,yc+k2/2);
  k4 = h*fname(tc+h,yc+k3);
  ynew = yc+(k1+2*k2+2*k3+k4)/6;
   
elseif k==5 then
  k1 = h*fc;
  k2 = h*fname(tc+h/4,yc+k1/4);
  k3 = h*fname(tc+3*h/8,yc+3/32*k1+9/32*k2);
  k4 = h*fname(tc+12/13*h,yc+1932/2197*k1-7200/2197*k2+7296/2197*k3);
  k5 = h*fname(tc+h,yc+439/216*k1-8*k2+3680/513*k3-845/4104*k4);
  k6 = h*fname(tc+1/2*h,yc-8/27*k1+2*k2-3544/2565*k3+1859/4104*k4-11/40*k5);
  ynew = yc+16/135*k1+6656/12825*k3+28561/56430*k4-9/50*k5+2/55*k6;
end
tnew = tc+h;
fnew = fname(tnew,ynew);