deff('y=f(t,x)','y(1)=-x(1)*x(2)/2+x(1)/2;y(2)=x(1)*x(2)-x(2)');
//
T=250;  // instant final
//
N=evstr(x_dialog('nombre de points N','100000'));
y0init=evstr(x_dialog('donnée initiale','[2;4]'));
//
dt=T/N;t=0:dt:T;
//
//
// solution Euler
//
y0=y0init
yeu=y0;
for i=1:N
  y0=y0+dt*f(t(i),y0);
  yeu=[yeu,y0]  
end
xset('window',2);
clf;
plot2d([t',t'],[yeu(1,:)',yeu(2,:)'],leg='proies@predateur')
xtitle('Euler')
//
xset('window',3);
clf
plot2d(yeu(1,:),yeu(2,:))
xtitle('Euler')
//
//
// solution RK4
//
//
y0=y0init
yRK=y0;
for i=1:N
//
t1=t(i);
t2=t(i)+dt/2;
t3=t(i)+dt/2;
t4=t(i)+dt;
//
y1=y0
y2=y0+dt*1/2*f(t1,y1);
y3=y0+dt*1/2*f(t2,y2);
y4=y0+dt*f(t3,y3);
//
y0=y0+dt*(1/6*f(t1,y1)+2/6*f(t2,y2)+2/6*f(t3,y3)+1/6*f(t4,y4))
//
  yRK=[yRK,y0]  
end
//
xset('window',0);
clf;
plot2d([t',t'],[yRK(1,:)',yRK(2,:)'],leg='proies@predateur')
xtitle('Runge Kutta 4')
//
xset('window',1);
clf
plot2d(yRK(1,:),yRK(2,:))
xtitle('Runge Kutta 4')
//

