function preda_predatore_rk4(x0, y0, t, dt, A, B, C, D, E)
% programma preda_predatore

% warning('off');

if nargin < 9, 
    disp(sprintf(' inserire x0, y0, t, dt, A, B, C, D, E'))
    return
end

P1=[0 0]';
P2=[B/D (A-E*B/D)/C]';
P3=[A/E 0]';
disp(sprintf('punti di equilibrio: P1=(%d,%d) P2=(%d,%d) P3=(%d,%d)',P1(1,1),P1(2,1),P2(1,1),P2(2,1),P3(1,1),P3(2,1)))
x=zeros(t/dt,1);
y=zeros(t/dt,1);
x(1)=x0;
y(1)=y0;

for i = 1:t/dt
    kx1=(A*x(i)-C*x(i)*y(i)-E*x(i)^2)*dt;
    ky1=(-B*y(i)+D*x(i)*y(i))*dt;
    xm=x(i)+kx1/2;
    ym=y(i)+ky1/2;
    kx2=(A*xm-C*xm*ym-E*xm^2)*dt;
    ky2=(-B*ym+D*xm*ym)*dt;
    xm=x(i)+kx2/2;
    ym=y(i)+ky2/2;
    kx3=(A*xm-C*xm*ym-E*xm^2)*dt;
    ky3=(-B*ym+D*xm*ym)*dt;
    xm=x(i)+kx3;
    ym=y(i)+ky3;
    kx4=(A*xm-C*xm*ym-E*xm^2)*dt;
    ky4=(-B*ym+D*xm*ym)*dt;
    x(i+1,1)=x(i)+kx1/6+kx2/3+kx3/3+kx4/6;
    y(i+1,1)=y(i)+ky1/6+ky2/3+ky3/3+ky4/6;
end


timef=size(x,1)*dt;
t=(0:dt:timef-dt)';

figure(1);clf;
 plot(x, y);
 hold on;
 plot(P1(1),P1(2), 'xr')
 hold on
 plot(P2(1),P2(2), '*r')
 hold on
 plot(P3(1),P3(2), '+r')
 hold on;
 plot(x0,y0,'og')
s_title=sprintf('Traiettoria popolazioni');
title(s_title);
ylabel('Predatori')
xlabel('Prede');
xlim([0 (max(x)+1)*1.2]);
ylim([0 (max(y)+1)*1.2]);


figure(2);clf;
  subplot(2,1,1)
  plot(t,x)
  ylabel('Popolazione prede')
  xlabel('tempo');
  xlim([t(1) t(end)]);
  ylim([0 (max(x)+1)*1.2]);
  s_title=sprintf('Prede-predatori');
  title(s_title);
  subplot(2,1,2)
  plot(t,y)
  ylabel('Popolazione predatori')
  xlabel('tempo');
  xlim([t(1) t(end)]);
  ylim([0 (max(y)+1)*1.2]);

return