clear all
%solving logistic with Euler
e=1;
a=0.1;
delta=0.01;
T=1;
Nt=ceil(T/delta);%floor
X=linspace(0,1.1,401);


A=linspace(0.03,0.3,101);
for ir=1:length(A)
    a=A(ir);


for ix=1:length(X)

    x0=X(ix);
x(1)=x0;
t(1)=0;
for it=2:Nt+1
    t(it)=(it-1)*delta;
    y=x(it-1);
    v=y.*(1-y)-a*(1+e*sin(2*pi*t(it-1)));
    x(it)=x(it-1)+delta*v;
end;

% figure(10)
% %plot(t,x,'b','linewidth',2,'markersize',16)
% plot(t,x,'--k','linewidth',2,'markersize',12)
% xlabel('$t$ ','FontSize',24,'Interpreter','latex');
% ylabel('$x(t)$ ','FontSize',24,'Interpreter','latex');
% set(gca,'Fontsize',24)
% grid on
% hold on
% drawnow
% pause(0.2)

   P(ix)=x(end);

end


% %figure(10);hold off
%  sd=min(abs(X-P))
%  %pause
%  if sd<1.e-3
%  X1=X(X<0.5);n1=length(X1);P1=P(1:n1);
%  X2=X(n1+1:end);P2=P(n1+1:end);
% [sd1,ind1]=min(abs(X1-P1));
% [sd2,ind2]=min(abs(X2-P2));
% xe1(ir)=X1(ind1);if xe1(ir)<0;xe1(ir)=NaN;end
% xe2(ir)=X2(ind2);if xe2(ir)<0;xe2(ir)=NaN;end
%  else
%    xe1(ir)=NaN;xe2(ir)=NaN;
%  end

[sd,ind]=min(abs(X-P));
sd=min(abs(X-P))
 if sd<1.e-3
X1=[X(1:ind-5),X(ind+5:end)];P1=[P(1:ind-5),P(ind+5:end)];
[sd1,ind1]=min(abs(X1-P1));
xe1(ir)=X(ind);if xe1(ir)<0;xe1(ir)=-5;end
xe2(ir)=X1(ind1);if xe2(ir)<0;xe2(ir)=-5;end
else
   xe1(ir)=NaN;xe2(ir)=NaN;
 end

figure(20)
%plot(t,x,'b','linewidth',2,'markersize',16)
plot(X,P,'b',X,X,'r',xe1(ir),xe1(ir),'or',xe2(ir),xe2(ir),'or','linewidth',2,'markersize',12)
xlabel('$x$ ','FontSize',24,'Interpreter','latex');
ylabel('$p(x)$ ','FontSize',24,'Interpreter','latex');
set(gca,'Fontsize',24)
grid on
drawnow
%hold on




end

% figure(13)
% plot(R,xe1,'--r',R,xe2,'r','linewidth',2,'markersize',16)
% xlabel('$r$ ','FontSize',24,'Interpreter','latex');
% ylabel('$x_e(r)$ ','FontSize',24,'Interpreter','latex');
% hold off

figure(130)
plot(A,xe1,'.r',A,xe2,'.r','linewidth',2,'markersize',16)
xlabel('$a$ ','FontSize',24,'Interpreter','latex');
ylabel('$x_e(a)$ ','FontSize',24,'Interpreter','latex');
hold off