function [h,t,u] = heun(fun, t0, y0, N, a)

h=...;
t=...;
u=zeros(length(y0),N+1); % List initialization

% Write your code here
% ...
end