clear all
global a e


a=3/8;

AA=linspace(0,1/4,100);
AA=1/8;

tic
for ic=1:length(AA);
    a=AA(ic);
e=1;

xp=1/2+sqrt(1/4-a);
xm=1/2-sqrt(1/4-a);


tf=20;
Nt=100;
T=linspace(0,tf,Nt);
dt=T(2)-T(1);
Xt=0*T;

nx=500;
X0=linspace(0,1.5,nx);
X0=0.2;

P=0*X0;
%nx=1;
%X0=0.51;


for ib=1:length(X0);
    
x0=X0(ib);

xn=x0;
Xt(1)=xn;

for ia=2:Nt
    
tn=T(ia-1);

k1=fish(xn,tn);
k2=fish(xn+dt*k1/2,tn+dt/2);
k3=fish(xn+dt*k2/2,tn+dt/2);
k4=fish(xn+dt*k3,tn+dt);
xn=xn+dt*(k1+2*k2+2*k3+k4)/6;

if xn <0
    xn=0;
end;

Xt(ia)=xn;

end;

P(ib)=Xt(ia);

figure(1)
plot(T,Xt,'r',T,0*T,'k',T,xp*exp(0*T),'--r',T,xm*exp(0*T),'--r','linewidth',2);
xlabel('t'),ylabel('x(t)');
drawnow
hold on

% figure(1)
% plot(T,Xt,'r','linewidth',2);
% xlabel('t'),ylabel('x(t)');
% drawnow
% hold on


end

% figure(2)
% plot(X0,P,'r',X0,X0,'--k','linewidth',2)
% xlabel('x');ylabel('p(x)')
% hold on

% 
% [sd,ind]=min(abs(P-X0));
% 
% xmin1=X0(ind);
% P1=P;
% P1(ind)=40;
% [sd,ind]=min(abs(P1-X0));
% 
% xmin2=X0(ind);
% 
% if xmin1<0.1;xmin1=0;end;
% if xmin2<0.1;xmin2=0;end;
% 
% B1(ic)=xmin1;
% B2(ic)=xmin2;
if abs(imag(xp))>0.00001;
    BC1(ic)=-1;
    BC2(ic)=-1;
else;
    BC1(ic)=xp;
    BC2(ic)=xm;
end;
iter=0;
for it=2:length(X0);
    
    gh=sign((X0(it)-P(it))*(X0(it-1)-P(it-1)));
    
    
    if gh <0 ;iter=iter+1;
        if iter==1;B1(ic)=(X0(it)+X0(it-1))/2;end;%X0(it-1)+P(it-1)*(X0(it)-X0(it-1))/(P(it)-P(it-1));end;
        if iter==2;B2(ic)=(X0(it)+X0(it-1))/2;end;%X0(it-1)+P(it-1)*(X0(it)-X0(it-1))/(P(it)-P(it-1));end;
    end;
    if iter==1; B2(ic)=-1;end;
    if iter==0;B1(ic)=-1;B2(ic)=-1;end;
        
end;

    
    




end;

toc
figure(3)
plot(AA,BC1,'r',AA,BC2,'b',AA,B1,'ok',AA,B2,'ok','linewidth',2)
axis([0 max(AA) 0 1.1])
xlabel('a');ylabel('x_e(a)')