function [ x, y ] = rk3 ( f_ode, xRange, yInitial, numSteps )
% [ x, y ] = rk2 ( f_ode, xRange, yInitial, numSteps )
% comments including the signature, meanings of variables,
% math methods, your name and the date

x(1,1) = xRange(1);
h = ( xRange(2) - xRange(1) ) / numSteps;
y(:,1) = yInitial;
for k = 1 : numSteps
  xa = x(k)+h/2 ;
  ya = y(:,k)+h/2*f_ode(x(k),y(:,k)) ;
  xb = x(k)+h;
  yb = y(:,k)+h*(2*f_ode(xa,ya)-f_ode(x(k),y(:,k)));
  x(1,k+1) = x(1,k) + h;
  y(:,k+1) = y(:,k) + h * (f_ode(x(k),y(:,k))+4*f_ode(xa,ya)+f_ode(xb,yb))/6;
end
