%% Heun method
% Solve the IVP y'(x)=f(x,y) with y(y0)=x0. The function f is given as a
% symbolic function.
% Given: step size h and number of steps. 

function [xH,yH] = Heun(f,x0,y0,h,N)

if length(h) == 1
    h = h*ones(1,N);
end
xH = zeros(1,N+1); yH = zeros(1,N+1);   % Initialization 
xH(1) = x0; yH(1) = y0;                 % Initial condition

for i=1:N                               % For loop over the steps
    k1 = f(xH(i),yH(i));                % k1
    k2 = f(xH(i)+h(i),yH(i)+h(i)*k1);         % k2
    k = (k1+k2)/2;                      % k
    xH(i+1) = xH(i)+h(i);                  % Next x
    yH(i+1) = yH(i)+h(i)*k;                % Heun iteration - next y
end