% 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