function y = ode_rk4 (f, int, y0, N)

  h = (int(2)-int(1))/N;
  y = zeros (length(y0), N+1);
  y(:,1) = y0;
  for j = 1:N
    tn = int(1) + j*h;
    yn = y(:,j);
    k1 = feval (f, tn, yn);
    k2 = feval (f, tn+h/2, yn + (h/2)*k1);
    k3 = feval (f, tn+h/2, yn + (h/2)*k2);
    k4 = feval (f, tn+h, yn + h*k3);
    y(:,j+1) = yn + (h/6)*(k1+2*k2+2*k3+k4);
  end

end