function y = RK44(x0, a, n, y0, f)
   h = a/n;
  y(1) = y0;
  x = x0 + h * [0:n];
  for i=1:n
    k1 = h*feval(f, x(i), y(i));
    k2=h*feval(f, x(i) + h/2, y(i) +k1/2);
    k3 = h*feval(f, x(i) + h/2, y(i) +k2/2);
    k4 = h*feval(f, x(i) + h, y(i) +k3);
    y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 +k4)/6;
  endfor
endfunction