function [xvals yvals]=RungeKuttaMethod(f,xmin,xmax,steps,y0)
  % calculate the step size
  h=(xmax-xmin)/steps;
  % now set up a list of x values
  % need the list to have steps+1 numbers in it
  % because the initial x is counted as 1
  xvals=linspace(xmin,xmax,steps+1);
  x=xmin;
  % initialize y and yvals
  y=y0;
  yvals=[y];
  % now a for loop to build the yvals
  for i=1:steps
      % x,y have current values
      k1=f(x,y);
      k2=f(x+h/2,y+h*k1/2);
      k3=f(x+h/2,y+h*k2/2);
      k4=f(x+h,y+h*k3);
      y=y+h/6*(k1+2*k2+2*k3+k4);
      % add y to the list of yvals
      yvals=[yvals y];
      % update x
      x=xvals(i+1);
  end
end
