//  Volterra avec 4 méthodes (Euler, heun, RK4 ou ode)
//
//
a=1;b=1;c=1;d=1;y01=3;y02=0.5
deff('ydot=fvol(t,y)','ydot=[a*y(1)-b*y(1)*y(2);-c*y(2)+d*y(1)*y(2)]')
//
N=evstr(x_dialog('nombre de points N','2000'));
Nm=evstr(x_dialog('choix de la méthode( 1: Euler, 2: Heun, 3: RK4, 4: ode)','1'));
//
T=40
h=T/N
t=0:h:T
//
y=[y01;y02]

if (Nm<4) then

Yapp=[y]
for i=1:N
//
if (Nm==1) then 
    y=y+h*fvol(t(i),y);// Euler
end
//
if (Nm==2) then     
y=y+1/2*h*(fvol(t(i),y)+fvol(t(i+1),y+h*fvol(t(i),y))); // Heun  
end
// 
if(Nm==3) then 

t1=t(i);
t2=t(i)+dt/2;
t3=t(i)+dt/2;
t4=t(i)+dt;
//
y1=y
y2=y+h*1/2*fvol(t1,y1);
y3=y+h*1/2*fvol(t2,y2);
y4=y+h*fvol(t3,y3);
//
y=y+h*(1/6*fvol(t1,y1)+2/6*fvol(t2,y2)+2/6*fvol(t3,y3)+1/6*fvol(t4,y4))
//    
end
//              
Yapp=[Yapp,y]
end
end
if (Nm==4)  then
  Yapp=ode(y,0,t,fvol); // solution ode
end  
//
figure(1)
clf()
plot(t,Yapp(1,:),'b-')
plot(t,Yapp(2,:),'black-')
//
figure(2)
clf()
plot(Yapp(1,:),Yapp(2,:))
