% Heun's method for numerical integration of ODEs, vector version
% INPUT:
% 1) an anonymous function f(t, x)
%    of two variables:
%    t: numerical
%    x: column vector
% 2) a vector T of t-positions (including t0 as the first one)
% 3) a column vector x0 of initial values of x for t=t0

% OUTPUT:
% a matrix of as many rows as x0 and as many columns as T above,
% which approximate the solution to the (vector) ODE given by f.
function [y] = heunvector(f, x, y0)
    % first of all, create 'y' with the adequate size
    y = zeros(length(y0),length(x));
    % and store the initial condition at the first place
    y(:,1) = y0;

    % Run Heun's loop
    for s = 2:length(x)
        % "velocity at this point"
        v1 = f(x(s-1), y(:,s-1));
        % "next point for Euler"  
        P = y(:,s-1) + (x(s) - x(s-1))*v1;
        % "velocity at this point"
        v2 = f(x(s), P);
        % mean velocity
        v  = (v1 + v2)/2;
        % next point
        y(:, s) = y(:,s-1) + (x(s) - x(s-1))*v;
    end
end