%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to solve the initial value problem
%
%   dy/dx = F(x,y), x = [a,b], y(a) = ya
%
% with the RK4 Runge-Kutta method. The vector u(2:N+1) approximates the true
% solution y(1:N) at x = a + ih, i = 1:N. Notice u(1) = y(a).
%
% With Matlab you need to define F with somename = @(x,y) expression before
% invoking this function (with arglist (...,somename)).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u = int_rk4(a, b, ya, N, F)

  h = (b-a) / N;
  u = zeros(1,N);

  u(1) = ya;

  for i = 1:N
    x = a + (i-1)*h;

    k1 = F(x, u(i));
    k2 = F(x+0.5*h, u(i)+0.5*h*k1);
    k3 = F(x+0.5*h, u(i)+0.5*h*k2);
    k4 = F(x+h, u(i)+h*k3);

    u(i+1) = u(i) + (h/6.)*(k1+2*k2+2*k3+k4);
end