//  comparaison méthode Euler et méthode de Heun, cas exact
//
C=1;D=4;y0=2;
//
deff('ydot=f(t,y)','ydot=C*y+D')
deff('y=yex(t)','y=-D/C+(y0+D/C)*exp(C*t)')
//
T=1;
errmaxe=[];errmaxh=[];
for j=1:5
N=10^j
h=T/N;
t=0:h:T;
//
ye=y0;
yh=y0;
//
Yerre=[0];Yerrh=[0];
for i=1:N
 ye=ye+h*f(t(i),ye); // Euler
 yh=yh+h/2*(f(t(i),yh)+f(t(i)+h,yh+h*f(t(i),yh))); // Heun
//
 Yerre=[Yerre,ye-yex(t(i+1))];
 Yerrh=[Yerrh,yh-yex(t(i+1))];
end
errmaxe=[errmaxe,max(abs(Yerre))];
errmaxh=[errmaxh,max(abs(Yerrh))];
//
end
clf()
plot2d(log(10.^(1:5)),log(errmaxe),1)
plot2d(log([10,100]),log([100,10]),2)
plot2d(log(10.^(1:5)),log(errmaxh),3)
plot2d(log([10,100]),log([100,1]),4)
legends(['erreur Euler';'ordre 1';'erreur Heun';'ordre 2'],[1,2,3,4])
