% syntbac
x1=0.0360;
x2=0.4764;
Tmax=100;
figure(1),
plot([0,1],[x2,x2],'r');
axis([-0.1,1.1,-0.1,1.1]);
hold on,
options = odeset('RelTol',1e-5,'AbsTol',[1e-4 1e-4 1e-4 1e-4],'Events',@OutOfSpace);
X0=[x1;x2;-1;0];
for k=1:100,
%     X0=[x1;x2;cos(2*pi*k/1000);sin(2*pi*k/1000)];
%     [T,X]=ode45(@optbac,[0,Tmax],X0,options);
%     X=real(X);
%     plot(X(:,1),X(:,2),'b');
    x1=(k/500)*(1-x2);
    X0=[x1;x2;0;1];
    [T,X]=ode45(@optbac,[0,Tmax],X0,options);
    X=real(X);
    plot(X(:,1),X(:,2),'b');
    if abs(X(end,1)-0.2809)<0.01,Xsaved=X0;end,
    X0=[x1;x2;0;-1];
    [T,X]=ode45(@optbac,[0,Tmax],X0,options);
    X=real(X);
    plot(X(:,1),X(:,2),'g');
    if abs(X(end,1)-0.2809)<0.01,Xsaved=X0;end,
end,

% x1=0.0360;
% x2=0.4764;
% Xsaved=[x1;x2;-1;0];
disp(Xsaved'),
X0=Xsaved;
[T,X]=ode45(@optbac,[0,Tmax],X0,options);
X=real(X);
x1=X(:,1);x2=X(:,2);x3=1-x1-x2;
p1=X(:,3);p2=X(:,4);
figure(1),hold on,
np=plot(x1,x2,'r');
set(np,'LineWidth',2);


hold off
zoom on,

figure(2),
subplot(2,1,1),
plot(T(end)-T(end:-1:1),x1(end:-1:1),'r',...
     T(end)-T(end:-1:1),x2(end:-1:1),'b',...
     T(end)-T(end:-1:1),x3(end:-1:1),'g'),title('Trajectories r/b/g');

subplot(2,2,3),plot(T(end)-T(end:-1:1),p1(end:-1:1),'r'),title('p1');

u=0.7*((p1<0).*(x1>=0));
H=p1.*(sqrt(x2)-u)+p2.*(sqrt(x3)-sqrt(x2));
subplot(2,2,4),plot(T(end)-T(end:-1:1),H(end:-1:1)),title('Hamiltonian');