% Heun's method for numerical integration of ODEs
% INPUT:
% 1) an anonymous function f
% 2) a vector of x-positions (including x0 as the first one)
% 3) an initial value y0, corresponding to x0

% OUTPUT:
% a vector y of values of the approximate solution at 'x'
function [y] = heun(f, x, y0)
    % first of all, create 'y' with the same length as x
    y = zeros(size(x));
    % and store the initial condition at the first place
    y(1) = y0;

    % Run Heun's loop
    for s = 2:length(x)
        k1 = f(x(s-1),y(s-1));
        y2 = y(s-1) + k1.*(x(s)-x(s-1));
        k2 = f(x(s), y2);
        k  = (k1+k2)./2;
        y(s) = y(s-1) + k*(x(s) - x(s-1));
    end
end