function dinamica_popolazione_lag(x0, A , B, tlag)
% programma dinamica_popolazione

if nargin < 4, 
    fprintf(' inserire x0, A, B, tlag')
    return
end


tmax=255;
DT=tmax/100;
N=floor(tmax/DT);
n=floor(tlag/DT);
y=zeros(N,1);
y(1)=x0;
for i = 2:n
     y(i)=x0;
end
y(n+1)=x0;
tic
for i = n:N
    y1=(A*y(i)-B*y(i)^2)*DT;
    yn=y(i)+y1/2;
    y2=(A*yn-B*yn^2)*DT;
    yn=y(i)+y2/2;
    y3=(A*yn-B*yn^2)*DT;
    yn=y(i)+y3;
    y4=(A*yn-B*yn^2)*DT;
    y(i+1,1)=y(i)+y1/6+y2/3+y3/3+y4/6;
end
toc

t=(0:DT:tmax);
figure(1);
plot(t, y);
ylabel('Popolazione')
xlabel('tempo');
s_title=sprintf('Dinamica popolazione RK4 lin.');
title(s_title);
xlim([t(1) t(end)]);
ylim([0 max(y)*1.2]);

figure(2);
semilogy(t, y);
ylabel('Popolazione')
xlabel('tempo');
s_title=sprintf('Dinamica popolazione RK4 log.');
title(s_title);
xlim([t(1) t(end)]);
ylim([0 max(y)*1.2]);

ylognorm=log(y/y(1));

t2=(0:DT*2:tmax);

yspline=spline(t,ylognorm,t2);

figure(3);
plot(t, ylognorm, 'b+', t2,yspline, 'r*');
ylabel('log(N/N0)')
xlabel('tempo');
s_title=sprintf('Dinamica popolazione RK4 log.');
title(s_title);
xlim([t(1) t(end)]);
ylim([0 max(ylognorm)*1.2]);


ymax=max(y)*1.7;
DY=ymax/500;
y=(0:DY:ymax)';
tic
U=-A/2*y.^2+B/3*y.^3;
toc
figure(4)
  plot(y,U)
  ylabel('Potenziale')
  xlabel('Popolazione');
  s_title=sprintf('Din. pop. potenziale efficace');
  title(s_title);

return