function y=rk_hw6(tmin, tmax, nsteps, y0)
%
%  User supplies initial time (tmin), final time (tmax),
%  number of time steps (nsteps), and vector of initial
%  values (y0).  (y0 must be a COLUMN vector.)
%  This code uses the fourth-order Runge-Kutta method
%  to solve the initial value problem
%
%               y'=f(t,y), y(tmin) = y0
% 
%  The code returns an array y of size (d,nsteps+1), where d
%  is the number of unknowns (the length of y0), containing the
%  values of y at each timestep (including the initial time).
%
%    Load initial values.
%
d = length(y0);
y = zeros(d,nsteps+1);
y(:,1) = y0;
yk = y0;
%
%    Loop over time steps.
%
h = (tmax-tmin)/nsteps;

for k=1:nsteps,
  tk = tmin + k*h;

  F1 = f(tk, yk);
  F2 = f(tk+h/2, yk+(h/2)*F1);
  F3 = f(tk+h/2, yk+(h/2)*F2);
  F4 = f(tk+h, yk+h*F3);

  yk = yk + (h/6)*(F1 + 2*F2 + 2*F3 + F4);
  y(:,k+1) = yk;
end;

%
%  Function defining y'.
%
function fty = f(t,y)
        fty(1,1) = -1000*y(1) + y(2);
        fty(2,1) = -0.1*y(2);
	return;
