function [x,y] = my_Heun2(fun,a,b,y0,N)
% Heun method
% num_ODEs_gr.pdf pages 27 & 28 (explicit method)
% Integration interval [a,b]
% Initial condition y0
% Number of (partition) points N+1

x = zeros(N+1,1) ;
y = zeros(N+1,1) ;

h = (b-a)/N ;
x(1) = a ;
y(1) = y0 ;
for n = 1:N
     x(n+1) = x(n) + h ;
     k1 = feval(fun,x(n),y(n)) ;
     k2 = feval(fun,x(n)+h, y(n) + h*k1) ;
     y(n+1) = y(n) + h*(k1+k2)/2 ;
end    

end

