//  Euler comparaison cas exact
//
C=3;D=4;y0=5;
//
deff('ydot=f(t,y)','ydot=C*y+D')
deff('y=yex(t)','y=-D/C+(y0+D/C)*exp(C*t)')
//
T=1;
Ntot=10^(2:0.5:5)
errtot_e=[];
errtot_h=[];
errtot_rk=[];
//
for n=1:length(Ntot)
 N=Ntot(n)
 h=T/N;
 t=0:h:T;
//
 ye=y0;
 yh=y0;
 yrk=y0;
//
 Yapp_e=[y0]
 Yapp_h=[y0]
 Yapp_rk=[y0]
//
 for i=1:N
 ye=ye+h*f(t(i),ye);
 yh=yh+h*(1/2*f(t(i),yh)+1/2*f(t(i)+h,yh+h*f(t(i),yh)))
//
t1=t(i);
t2=t(i)+h/2;
t3=t(i)+h/2;
t4=t(i)+h;
//
y1=yrk
y2=yrk+h*1/2*f(t1,y1);
y3=yrk+h*1/2*f(t2,y2);
y4=yrk+h*f(t3,y3);
//
yrk=yrk+h*(1/6*f(t1,y1)+2/6*f(t2,y2)+2/6*f(t3,y3)+1/6*f(t4,y4))
//
Yapp_e=[Yapp_e,ye]
Yapp_h=[Yapp_h,yh]
Yapp_rk=[Yapp_rk,yrk]
 end
//figure(1)
//clf()
//plot(t,Yapp_e,'b-*');
//
 Yexacte=yex(t);
// plot(t,Yexacte,'black-')
erreur_e=max(abs(Yapp_e-Yexacte))
erreur_h=max(abs(Yapp_h-Yexacte))
erreur_rk=max(abs(Yapp_rk-Yexacte))
//
//disp(erreur_e)
//disp(erreur_h)
//
errtot_e=[errtot_e,erreur_e]
errtot_h=[errtot_h,erreur_h]
errtot_rk=[errtot_rk,erreur_rk]
end
clf()
plot(1/log(10)*log(Ntot),1/log(10)*log(errtot_e),'b-*')
plot(1/log(10)*log(Ntot),1/log(10)*log(errtot_h),'black*-')
plot(1/log(10)*log(Ntot),1/log(10)*log(errtot_rk),'red*-')
plot(1/log(10)*log([100,1000]),1/log(10)*log([1/100,1/1000]),'y-')
plot(1/log(10)*log([100,1000]),1/log(10)*log([1/100,1/10000]),'y-')
plot(1/log(10)*log([100,1000]),1/log(10)*log([1/100,1/1000000]),'y-')
