function [t,u]=Heun(f,tI,tF,yI,N)
h=(tF-tI)/N;
t=[tI:h:tF];
u(:,1)=yI;
for n=1:N
    K1=feval(f,t(n),u(:,n));
    K2=feval(f,t(n)+h,u(:,n)+h*K1);
    u(:,n+1)=u(:,n)+h/2*(K1+K2);
end