clear all

global h e

h=1/8; e=0.5;


H=linspace(0.02,1/4+0.2,201);
%H=[1/8,1/2]
Nt=101;
tt=linspace(0,2*pi,Nt);
dt=tt(2)-tt(1);

X0=linspace(0,2,301);
P=0*X0;
XF1=0*H;XF2=0*H;


for ih=1:length(H)
   h=H(ih);

iter=0;
iflag=0;
for ia=1:length(X0)

x0=X0(ia);

xt(1)=x0;

for it=2:Nt
    ta=tt(it);
    tprev=tt(it-1);
    xprev=xt(it-1);

    k1=dt*fishdot(xprev,tprev);
    k2=dt*fishdot(xprev+k1/2,tprev+dt/2);
    k3=dt*fishdot(xprev+k2/2,tprev+dt/2);
    k4=dt*fishdot(xprev+k3,tprev+dt);

    xt(it)=xprev+(1/6)*(k1+2*k2+2*k3+k4);
    if xt(it)<0
        xt(it)=0;
    end

end

% figure(400)
% plot(tt,xt,'b','linewidth',2)
% xlabel('$t$','interpreter','latex')
% ylabel('$x(t)$','interpreter','latex')
% set(gca,'fontsize',18)
% hold on
% drawnow


P(ia)=xt(it);

f(ia)=P(ia)-x0;




if ia>1
g=f(ia)*f(ia-1);



if iter==0
    if g<0
        xf1=x0;
        iter=iter+1;
    end
end

if iter==1
    if g<0
        xf2=x0;
        iflag=1;
    end
end


end
if iflag==0
    xf1=NaN;xf2=NaN;
end;


end

figure(501);
plot(X0,P,'r',X0,X0,'k',xf1,xf1,'ob',xf2,xf2,'ob','linewidth',2,'markersize',15)
xlabel('$x$','interpreter','latex')
ylabel('$P(x)$','interpreter','latex')
set(gca,'fontsize',18)
drawnow

XF1(ih)=xf1;
XF2(ih)=xf2;

end

figure(50);hold on
plot(H,XF1,'--k',H,XF2,'k','linewidth',2,'markersize',15)
xlabel('$h$','interpreter','latex')
ylabel('$x_{eq}$','interpreter','latex')
xlim([0,H(end)])
set(gca,'fontsize',18)
drawnow









