function lorentz(x0, y0, z0, t, dt, B, C, D)
% programma preda_predatore

% warning('off');

if nargin < 8, 
    fprintf(' inserire x0, y0, z0, t, dt, B, C, D\n')
    return
end

P1=[0 0 0]';
P2=[sqrt(B*(D-1)) sqrt(B*(D-1)) D-1]';
P3=[-sqrt(B*(D-1)) -sqrt(B*(D-1)) D-1]';
fprintf('punti di equilibrio:\n P1=(%.3f,%.3f,%.3f)\n P2=(%.3f,%.3f,%.3f)\n P3=(%.3f,%.3f,%.3f)\n',P1(1,1),P1(2,1),P1(3,1),P2(1,1),P2(2,1),P2(3,1),P3(1,1),P3(2,1),P3(3,1))
x=zeros(t/dt,1);
y=zeros(t/dt,1);
z=zeros(t/dt,1);
x(1)=x0;
y(1)=y0;
z(1)=z0;

for i = 1:t/dt
    kx1=(-C*x(i)+C*y(i))*dt;
    ky1=(D*x(i)-y(i)-x(i)*z(i))*dt;
    kz1=(-B*z(i)+x(i)*y(i))*dt;
    xm=x(i)+kx1/2;
    ym=y(i)+ky1/2;
    zm=z(i)+kz1/2;
    kx2=(-C*xm+C*ym)*dt;
    ky2=(D*xm-ym-xm*zm)*dt;
    kz2=(-B*zm+xm*ym)*dt;
    xm=x(i)+kx2/2;
    ym=y(i)+ky2/2;
    zm=z(i)+kz2/2;
    kx3=(-C*xm+C*ym)*dt;
    ky3=(D*xm-ym-xm*zm)*dt;
    kz3=(-B*zm+xm*ym)*dt;
    xm=x(i)+kx3;
    ym=y(i)+ky3;
    zm=z(i)+kz3/2;
    kx4=(-C*xm+C*ym)*dt;
    ky4=(D*xm-ym-xm*zm)*dt;
    kz4=(-B*zm+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;
    z(i+1,1)=z(i)+kz1/6+kz2/3+kz3/3+kz4/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 lorentz xy');
title(s_title);
ylabel('y')
xlabel('x');
%xlim([0 (max(x)+1)*1.2]);
%ylim([0 (max(y)+1)*1.2]);


figure(2);clf;
  subplot(3,1,1)
  plot(t,x)
  ylabel('x')
  xlabel('tempo');
  xlim([t(1) t(end)]);
  %ylim([0 (max(x)+1)*1.2]);
  s_title=sprintf('velocità x y z');
  title(s_title);
  subplot(3,1,2)
  plot(t,y)
  ylabel('y')
  xlabel('tempo');
  xlim([t(1) t(end)]);
  %ylim([0 (max(y)+1)*1.2]);
  subplot(3,1,3)
  plot(t,z)
  ylabel('z')
  xlabel('tempo');
  xlim([t(1) t(end)]);
  %ylim([0 (max(z)+1)*1.2]);

return