%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to solve the initial value problem
%
%   dy/dx = F(x,y), x = [a,b], y(a) = ya
%
% with the RK2 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_rk2(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+h, u(i)+h*k1);

    u(i+1) = u(i) + 0.5*h*(k1+k2);
end