%-------------------------------------------------------------------------%
%         This software is licensed by Creative Commons BY-NC-SA:         %
%      http://creativecommons.org/licenses/by-nc-sa/3.0/it/legalcode      %
%-------------------------------------------------------------------------%
%
% File: SeparatrixPoints(p1,tmax,tol,npi,par,E,nsp)
%
% Goal: script that finds a separatrix point
%
% Inputs: p1:    matrix 3x2 containig two initial conditions
%         tmax:  integration time
%         tol:   tolerance used during the bisection algorithm
%         npi:   number of equilibrium points to be interpolated
%         par:   model parameters
%         E:     matrix of the equilibria Mx3, where M is the number of
%                equilibria. The matrix is organized as follows:
%                row 1: a stable equilibrium point
%                from row 2 to row npi+1: equilibrium points to be 
%                interpolated
%                row npi+2: a stable equilibrium point
%                row npi+3: a stable equilibrium point
%                from row npi+4 to row M: remaining equilibria
%         nsp:   number of stable equilibrium points
%
% Outputs: point1, point2, point3: a point lying on the a separatrix
%                                  surface
%
%-------------------------------------------------------------------------%
function  [point1 point2 point3] = SeparatrixPoints(p1,tmax,tol,npi,...
    par,E,nsp)
% Declare options used in the integration routine
comps = [1,3]; options = odeset('NonNegative',comps);
% Integrate the system with the two initial conditions
tspan = [0,tmax];
y0 = p1(:,1); [T,Z1] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
y0 = p1(:,2); [T,Z2] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
% Fix the tolerances used to check where the initial conditions stabilize
size_E = size(E,1); tol1 = ones(size_E); tol2 = ones(size_E);
tol3 = ones(size_E); tol4 = ones(size_E);
for i = 1:size_E
    if (E(i,1) < 0.01)
        tol1(i) = 0.01;
    else
        tol1(i) = 10^(-1)*abs(E(i,1));
    end
end
for i = 1:size_E
    if (E(i,2) < 0.01)
        tol2(i) = 0.01;
    else
        tol2(i) = 10^(-1)*abs(E(i,2));
    end
end
for i = 1:size_E
    if (E(i,3) < 0.01)
        tol3(i) = 0.01;
    else
        tol3(i) = 10^(-1)*abs(E(i,3));
    end
end
% Fix the tolerances used to check if the two initial conditions stabilize
% or not towards different stable equilibria
tol0(1) = max(tol1(1),tol1(npi+2));
tol0(2) = max(tol2(1),tol2(npi+2));
tol0(3) = max(tol3(1),tol3(npi+2));
tol4(1) = max(tol0(1),tol1(npi+3));
tol4(2) = max(tol0(2),tol2(npi+3));
tol4(3) = max(tol0(3),tol3(npi+3));
tol4a(1) = max(tol1(1),tol1(npi+2));
tol4a(2) = max(tol2(1),tol2(npi+2));
tol4a(3) = max(tol3(1),tol3(npi+2));
tol4b(1) = max(tol1(1),tol1(npi+3));
tol4b(2) = max(tol2(1),tol2(npi+3));
tol4b(3) = max(tol3(1),tol3(npi+3));
% Check if the integration time is sufficient for both initial conditions
q = 1; flag = zeros(1,size_E);
for i = 1:size_E
    if ((abs(Z1(end,1) - E(i,1)) > tol1(i)) || ...
            (abs(Z1(end,2) - E(i,2)) > tol2(i)) || ...
            (abs(Z1(end,3) - E(i,3)) > tol3(i)))
        t(q) = 0; q = q + 1;
    else
        t(q) = 1; q = q + 1;
    end
end
% If the integrarion time is not sufficient return
if t == flag
    point1 = [-1 -1 -1]; point2 = [-1 -1 -1]; point3 = [-1 -1 -1];
    return
end
q = 1; flag = zeros(1,size_E);
for i = 1:size_E
    if ((abs(Z2(end,1) - E(i,1)) > tol1(i)) || ...
            (abs(Z2(end,2) - E(i,2)) > tol2(i)) || ...
            (abs(Z2(end,3) - E(i,3)) > tol3(i)))
        t(q) = 0; q = q + 1;
    else
        t(q) = 1; q = q + 1;
    end
end
% If the integrarion time is not sufficient return
if t == flag
    point1 = [-1 -1 -1]; point2 = [-1 -1 -1]; point3 = [-1 -1 -1];
    return
end
% Move the first initial condition if lies on a stable manifold of a
% saddle point or coincides with an equilibrium point
if ((abs(Z1(end,1) - E(1,1)) > tol1(1)) || ...
        (abs(Z1(end,2) - E(1,2)) > tol2(1)) || ...
        (abs(Z1(end,3) - E(1,3)) > tol3(1))) && ...
        ((abs(Z1(end,1) - E(npi+2,1)) > tol1(npi+2)) || ...
        (abs(Z1(end,2) - E(npi+2,2)) > tol2(npi+2)) || ...
        (abs(Z1(end,3) - E(npi+2,3)) > tol3(npi+2))) && ...
        ((abs(Z1(end,1) - E(npi+3,1)) > tol1(npi+3)) || ...
        (abs(Z1(end,2) - E(npi+3,2)) > tol2(npi+3)) || ...
        (abs(Z1(end,3) - E(npi+3,3)) > tol3(npi+3)))
    p2(1,1) = abs(p1(1,1) - (tol/4)); p2(2,1) = abs(p1(2,1));
    p2(3,1) = abs(p1(3,1));
    p3(1,1) = abs(p1(1,1)); p3(2,1) = abs(p1(2,1) - (tol/4));
    p3(3,1) = abs(p1(3,1));
    p4(1,1) = abs(p1(1,1)); p4(2,1) = abs(p1(2,1));
    p4(3,1) = abs(p1(3,1) - (tol/4));
    % Integrate the system with the new initial conditions
    y0 = p2(:,1); [T,Z1] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
    y0 = p3(:,1); [T,Z3] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
    y0 = p4(:,1); [T,Z4] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
    % Check if with the new initial conditions the integration time is
    % sufficient
    q = 1; flag = zeros(1,size_E);
    for i = 1:size_E
        if  ((abs(Z1(end,1) - E(i,1)) > tol1(i)) || ...
                (abs(Z1(end,2) - E(i,2)) > tol2(i)) || ...
                (abs(Z1(end,3) - E(i,3)) > tol3(i)))
            t(q) = 0; q = q + 1;
        else
            t(q) = 1; q = q + 1;
        end
    end
    if t == flag
        point1 = [-1 -1 -1]; point2 = [-1 -1 -1]; point3 = [-1 -1 -1];
        return
    end
    q = 1; flag = zeros(1,size_E);
    for i = 1:size_E
        if  ((abs(Z3(end,1) - E(i,1)) > tol1(i)) || ...
                (abs(Z3(end,2) - E(i,2)) > tol2(i))|| ...
                (abs(Z3(end,3) - E(i,3)) > tol3(i)))
            t(q) = 0; q = q + 1;
        else
            t(q) = 1; q = q + 1;
        end
    end
    if t == flag
        point1 = [-1 -1 -1]; point2 = [-1 -1 -1]; point3 = [-1 -1 -1];
        return
    end
    q = 1; flag = zeros(1,size_E);
    for i = 1:size_E
        if  ((abs(Z4(end,1) - E(i,1)) > tol1(i)) || ...
                (abs(Z4(end,2) - E(i,2)) > tol2(i)) || ...
                (abs(Z4(end,3) - E(i,3)) > tol3(i)))
            t(q) = 0; q = q + 1;
        else
            t(q) = 1; q = q + 1;
        end
    end
    if t == flag
        point1 = [-1 -1 -1]; point2 = [-1 -1 -1]; point3 = [-1 -1 -1];
        return
    end
    % Choose a new initial condition in the way that the two initial
    % conditions evolve towards different equilibria
    if (((abs(Z1(end,1) - E(1,1)) < tol1(1)) && ...
            (abs(Z1(end,2) - E(1,2)) < tol2(1)) &&...
            (abs(Z1(end,3) - E(1,3)) < tol3(1))) || ...
            ((abs(Z1(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
            (abs(Z1(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
            (abs(Z1(end,3) - E(npi+2,3)) < tol3(npi+2))) || ...
            ((abs(Z1(end,1) - E(npi+3,1)) < tol1(npi+3)) && ...
            (abs(Z1(end,2) - E(npi+2,3)) < tol2(npi+3)) && ...
            (abs(Z1(end,3) - E(npi+3,3)) < tol3(npi+3)))) &&...
            ((abs(Z1(end,1) - Z2(end,1)) > tol4(1))  || ...
            (abs(Z1(end,2) - Z2(end,2)) > tol4(2)) || ...
            (abs(Z1(end,3) - Z2(end,3)) > tol4(3)))
        p1(:,1) = p2(:,1);
    elseif  (((abs(Z3(end,1) - E(1,1)) < tol1(1)) && ...
            (abs(Z3(end,2) - E(1,2)) < tol2(1)) && ...
            (abs(Z3(end,3) - E(1,3)) < tol3(1))) || ...
            ((abs(Z3(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
            (abs(Z3(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
            (abs(Z3(end,3) - E(npi+2,3)) < tol3(npi+2))) ||...
            ((abs(Z3(end,1) - E(npi+3,1)) < tol1(npi+3)) && ...
            (abs(Z3(end,2) - E(npi+2,3)) < tol2(npi+3)) && ...
            (abs(Z3(end,3) - E(npi+3,3)) < tol3(npi+3)))) &&...
            ((abs(Z3(end,1) - Z2(end,1)) > tol4(1)) || ...
            (abs(Z3(end,2) - Z2(end,2)) > tol4(2)) || ...
            (abs(Z3(end,3) - Z2(end,3)) > tol4(3)))
        p1(:,1) = p3(:,1);
    elseif  (((abs(Z4(end,1) - E(1,1)) < tol1(1)) && ...
            (abs(Z4(end,2) - E(1,2)) < tol2(1)) && ...
            (abs(Z4(end,3) - E(1,3)) < tol3(1))) || ...
            ((abs(Z4(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
            (abs(Z4(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
            (abs(Z4(end,3) - E(npi+2,3)) < tol3(npi+2))) ||...
            ((abs(Z4(end,1) - E(npi+3,1)) < tol1(npi+3)) && ...
            (abs(Z4(end,2) - E(npi+2,3)) < tol2(npi+3)) && ...
            (abs(Z4(end,3) - E(npi+3,3)) < tol3(npi+3)))) && ...
            ((abs(Z4(end,1) - Z2(end,1)) > tol4(1))  || ...
            (abs(Z4(end,2) - Z2(end,2)) > tol4(2)) || ...
            (abs(Z4(end,3) - Z2(end,3)) > tol4(3)))
        p1(:,1) = p4(:,1);
    else
        % For the given initial conditions is not possible to performe 
        % bisection
        point1 = [0 0 0]; point2 = [0 0 0]; point3 = [0 0 0];
        return
    end
end
% Move the second initial condition if lies on a stable manifold of a
% saddle point or coincides with an equilibrium point
y0 = p1(:,1);
[T,Z1] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
if ((abs(Z2(end,1) - E(1,1)) > tol1(1)) || ...
        (abs(Z2(end,2) - E(1,2)) > tol2(1)) || ...
        (abs(Z2(end,3) - E(1,3)) > tol3(1))) && ...
        ((abs(Z2(end,1) - E(npi+2,1)) > tol1(npi+2)) || ...
        (abs(Z2(end,2) - E(npi+2,2)) > tol2(npi+2)) || ...
        (abs(Z2(end,3) - E(npi+2,3)) > tol3(npi+2))) && ...
        ((abs(Z2(end,1) - E(npi+3,1)) > tol1(npi+3)) ||...
        (abs(Z2(end,2) - E(npi+3,2)) > tol2(npi+3)) || ...
        (abs(Z2(end,3) - E(npi+3,3)) > tol3(npi+3)))
    p2(1,1) = abs(p1(1,2) - (tol/4)); p2(2,1) = abs(p1(2,2));
    p2(3,1) = abs(p1(3,2));
    p3(1,1) = abs(p1(1,2)); p3(2,1) = abs(p1(2,2) - (tol/4));
    p3(3,1) = abs(p1(3,2));
    p4(1,1) = abs(p1(1,2)); p4(2,1) = abs(p1(2,2));
    p4(3,1) = abs(p1(3,2) - (tol/4));
    % Integrate the system with the new initial conditions
    y0 = p2(:,1); [T,Z2] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
    y0 = p3(:,1); [T,Z3] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
    y0 = p4(:,1); [T,Z4] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
    % Check if with the new initial conditions the integration  time is
    % sufficient
    q = 1;
    flag = zeros(1,size_E);
    for i = 1:size_E
        if  ((abs(Z2(end,1) - E(i,1)) > tol1(i)) || ...
                (abs(Z2(end,2) - E(i,2)) > tol2(i)) || ...
                (abs(Z2(end,3) - E(i,3)) > tol3(i)))
            t(q) = 0; q = q + 1;
        else
            t(q) = 1; q = q + 1;
        end
    end
    if t == flag
        point1 = [-1 -1 -1]; point2 = [-1 -1 -1]; point3 = [-1 -1 -1];
        return
    end
    q = 1; flag = zeros(1,size_E);
    for i = 1:size_E
        if  ((abs(Z3(end,1) - E(i,1)) > tol1(i)) || ...
                (abs(Z3(end,2) - E(i,2)) > tol2(i)) || ...
                (abs(Z3(end,3) - E(i,3)) > tol3(i)))
            t(q) = 0; q = q + 1;
        else
            t(q) = 1; q = q + 1;
        end
    end
    if t == flag
        point1 = [-1 -1 -1]; point2 = [-1 -1 -1]; point3 = [-1 -1 -1];
        return
    end
    q = 1; flag = zeros(1,size_E);
    for i = 1:size_E
        if  ((abs(Z4(end,1) - E(i,1)) > tol1(i)) || ...
                (abs(Z4(end,2) - E(i,2)) > tol2(i)) || ...
                (abs(Z4(end,3) - E(i,3)) > tol3(i)))
            t(q) = 0; q = q + 1;
        else
            t(q) = 1; q = q + 1;
        end
    end
    if t == flag
        point1 = [-1 -1 -1]; point2 = [-1 -1 -1]; point3 = [-1 -1 -1];
        return
    end
    % Choose a new initial condition in the way that the two initial
    % conditions evolve towards different stable equilibria
    if (((abs(Z2(end,1) - E(1,1)) < tol1(1)) && ...
            (abs(Z2(end,2) - E(1,2)) < tol2(1)) && ...
            (abs(Z2(end,3) - E(1,3)) < tol3(1))) || ...
            ((abs(Z2(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
            (abs(Z2(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
            (abs(Z2(end,3) - E(npi+2,3)) < tol3(npi+2))) || ...
            ((abs(Z2(end,1) - E(npi+3,1)) < tol1(npi+3)) && ...
            (abs(Z2(end,2) - E(npi+2,3)) < tol2(npi+3)) && ...
            (abs(Z2(end,3) - E(npi+3,3)) < tol3(npi+3)))) &&...
            ((abs(Z1(end,1) - Z2(end,1)) > tol4(1))  || ...
            (abs(Z1(end,2) - Z2(end,2)) > tol4(2)) || ...
            (abs(Z1(end,3) - Z2(end,3)) > tol4(3)))
        p1(:,2) = p2(:,1);
    elseif  (((abs(Z3(end,1) - E(1,1)) < tol1(1)) && ...
            (abs(Z3(end,2) - E(1,2)) < tol2(1)) && ...
            (abs(Z3(end,3) - E(1,3)) < tol3(1))) || ...
            ((abs(Z3(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
            (abs(Z3(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
            (abs(Z3(end,3) - E(npi+2,3)) < tol3(npi+2))) || ...
            ((abs(Z3(end,1) - E(npi+3,1)) < tol1(npi+3)) && ...
            (abs(Z3(end,2) - E(npi+2,3)) < tol2(npi+3)) && ...
            (abs(Z3(end,3) - E(npi+3,3)) < tol3(npi+3)))) && ...
            ((abs(Z3(end,1) - Z1(end,1)) > tol4(1))  || ...
            (abs(Z3(end,2) - Z1(end,2)) > tol4(2)) || ...
            (abs(Z3(end,3) - Z1(end,3)) > tol4(3)))
        p1(:,2) = p3(:,1);
    elseif (((abs(Z4(end,1) - E(1,1)) < tol1(1)) && ...
            (abs(Z4(end,2) - E(1,2)) < tol2(1)) && ...
            (abs(Z4(end,3) - E(1,3)) < tol3(1))) || ...
            ((abs(Z4(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
            (abs(Z4(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
            (abs(Z4(end,3) - E(npi+2,3)) < tol3(npi+2))) || ...
            ((abs(Z4(end,1) - E(npi+3,1)) < tol1(npi+3)) && ...
            (abs(Z4(end,2) - E(npi+2,3)) < tol2(npi+3)) && ...
            (abs(Z4(end,3) - E(npi+3,3)) < tol3(npi+3)))) && ...
            ((abs(Z4(end,1) - Z1(end,1)) > tol4(1))  || ...
            (abs(Z4(end,2) - Z1(end,2)) > tol4(2)) || ...
            (abs(Z4(end,3) - Z1(end,3)) > tol4(3)))
        p1(:,2) = p4(:,1);
    else
        % For the given initial conditions is not possible to performe 
        % bisection
        point1 = [0 0 0]; point2 = [0 0 0]; point3 = [0 0 0];
        return
    end
end
y0 = p1(:,2); [T,Z2] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
% Check if one of the two initial condition evolves towards the first
% stable equilibrium point
if ((abs(Z1(end,1) - E(1,1)) < tol1(1)) && ...
        (abs(Z1(end,2) - E(1,2)) < tol2(1)) && ...
        (abs(Z1(end,3) - E(1,3)) < tol3(1))) || ...
        ((abs(Z2(end,1) - E(1,1)) < tol1(1)) && ...
        (abs(Z2(end,2) - E(1,2)) < tol2(1)) && ...
        (abs(Z2(end,3) - E(1,3)) < tol3(1)))
    % Check if the first initial condition evolves towards the first
    % stable equilibrium point
    if (abs(Z1(end,1) - E(1,1)) < tol1(1)) && ...
            (abs(Z1(end,2) - E(1,2)) < tol2(1)) && ...
            (abs(Z1(end,3) - E(1,3)) < tol3(1))
        p0 = p1;
    else
        % Rearrange the matrix of initial conditions such that the first
        % initial condition evolves towards the first stable equilibrium
        p0(1,1) = p1(1,2); p0(1,2) = p1(1,1); p0(2,1) = p1(2,2); 
        p0(2,2) = p1(2,1); p0(3,1) = p1(3,2); p0(3,2) = p1(3,1);
        p0 = [p0(1,1) p0(1,2); p0(2,1) p0(2,2); p0(3,1) p0(3,2)];
    end
    % Integrate the system with the new initial conditions
    y0 = p0(:,1); [T,Z1] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
    y0 = p0(:,2); [T,Z2] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
    % Check if the initial conditions evolve towards two different
    % stable equilibrium points
    if (abs(Z1(end,1) - Z2(end,1)) < tol4a(1)) && ...
            (abs(Z1(end,2) - Z2(end,2)) < tol4a(2)) && ...
            (abs(Z1(end,3) - Z2(end,3)) < tol4a(3))
        point1 = [0 0 0]; point2 = [0 0 0]; point3 = [0 0 0];
        return
    end
    % Check if the second initial condition evolves towards the second
    % stable equilibrium point
    if (abs(Z2(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
            (abs(Z2(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
            (abs(Z2(end,3) - E(npi+2,3)) < tol3(npi+2))
        % Initialize the matrix of trajectories
        D = [Z1(end,:);Z2(end,:)];
        % Initialize the matrix of mid points
        M = [p0(:,1)';p0(:,2)'];
        % Compute the first mid point
        mx = (p0(1,1) + p0(1,2))/2; my = (p0(2,1) + p0(2,2))/2;
        mz = (p0(3,1) + p0(3,2))/2;
        m = [mx my mz]';
        % Integrate the system with the mid point as initial condition
        [T,W] = ode45(@ProgIntegr3D,tspan,m,[],par,options);
        flag1 = 0;
        if nsp == 3
            % If the mid point lies in the basin of attraction of 
            % the third equilibrium point the routine stops
            if (abs(W(end,1) - E(npi+3,1)) < tol1(npi+3)) ...
                    && (abs(W(end,2) - E(npi+3,2)) < tol2(npi+3)) ...
                    && (abs(W(end,3) - E(npi+3,3)) < tol3(npi+3))
                point1 = [0 0 0]; point2 = [0 0 0]; point3 = [0 0 0];
                return
            end
        end
        % Update matrices
        D = [D; W(end,:)]; M = [M;m'];
        % Calculate the distance between the two initial conditions
        dist = sqrt((p0(1,1) - p0(1,2))^2 + (p0(2,1) - p0(2,2))^2 + ...
            (p0(3,1) - p0(3,2))^2);
        % Bisection routine
        while (dist > tol)
            if (abs(W(end,1) - E(1,1)) < tol1(1)) && ...
                    (abs(W(end,2) - E(1,2)) < tol2(1)) && ...
                    (abs(W(end,3) - E(1,3)) < tol3(1))
                p0(:,1) = m;
            else
                p0(:,2) = m;
            end
            mx = (p0(1,1) + p0(1,2))/2; my = (p0(2,1) + p0(2,2))/2;
            mz = (p0(3,1) + p0(3,2))/2;
            m = [mx my mz]';
            [T,W] = ode45(@ProgIntegr3D,tspan,m,[],par,options);
            flag1 = 0;
            if nsp == 3
                % If the mid point lies in the basin of attraction of
                % the third equilibrium point the routine stops
                if (abs(W(end,1) - E(npi+3,1)) < tol1(npi+3)) && ...
                        (abs(W(end,2) - E(npi+3,2)) < tol2(npi+3)) && ...
                        (abs(W(end,3) - E(npi+3,3)) < tol3(npi+3))
                    flag1 = 1;
                    break
                end
            end
            dist =(sqrt((p0(1,1) - p0(1,2))^2 + (p0(2,1) - p0(2,2))^2 +  ...
                (p0(3,1) - p0(3,2))^2));
            b = W(end,:); D = [D;b]; M = [M;m'];
        end
        if flag1 == 1
            point1 = [0 0 0]; point2 = [0 0 0]; point3 = [0 0 0];
            return
        end
        j1 = 0; j = size(D,1);
        while (j1 == 0)
            % Find in D the index of last trajectory evolving towards the
            % first stable equilibrium
            if (abs(D(j,1) - E(1,1)) < tol1(1)) && ...
                    (abs(D(j,2) - E(1,2)) < tol2(1)) && ...
                    (abs(D(j,3) - E(1,3)) < tol3(1))
                r = j; j1 = 1;
            else
                j = j - 1;
            end
        end
        % Find in D the index of last trajectory evolving towards the
        % second stable equilibrium
        j1 = 0; j = size(D,1);
        while (j1 == 0)
            if (abs(D(j,1) - E(1,1)) > tol1(1)) || ...
                    (abs(D(j,2) - E(1,2)) > tol2(1)) || ...
                    (abs(D(j,3) - E(1,3)) > tol3(1))
                s = j; j1 = 1;
            else
                j = j - 1;
            end
        end
        % Declare the separatrix point
        point1(1) = (M(r,1) + M(s,1))/2; point1(2) = (M(r,2) + M(s,2))/2;
        point1(3) = (M(r,3) + M(s,3))/2;
        point2(1) = 0; point2(2) = 0; point2(3) = 0;
        point3(1) = 0; point3(2) = 0; point3(3) = 0;
        return
    else
        % The second initial condition evolves towards the third
        % stable equilibrium
        % Initialize the matrix of trajectories
        D = [Z1(end,:);Z2(end,:)];
        % Initialize the matrix of mid points
        M = [p0(:,1)';p0(:,2)'];
        % Compute the first mid point
        mx = (p0(1,1) + p0(1,2))/2; my = (p0(2,1) + p0(2,2))/2;
        mz = (p0(3,1) + p0(3,2))/2;
        m = [mx my mz]';
        % Integrate the system with the mid point as initial condition
        [T,W] = ode45(@ProgIntegr3D,tspan,m,[],par,options);
        flag1 = 0;
        if nsp == 3
            % If the mid point lies in the basin of attraction of the
            % second stable equilibrium point the routine stops
            if (abs(W(end,1) - E(npi+2,1)) < tol1(npi+2)) ...
                    && (abs(W(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
                    (abs(W(end,3) - E(npi+2,3)) < tol3(npi+2))
                point1 = [0 0 0]; point2 = [0 0 0]; point3 = [0 0 0];
                return
            end
        end
        % Update matrices
        D = [D; W(end,:)]; M = [M;m'];
        dist = sqrt((p0(1,1) - p0(1,2))^2 + (p0(2,1) - p0(2,2))^2 + ...
            (p0(3,1) - p0(3,2))^2);
        % Bisection routine
        while (dist > tol)
            if (abs(W(end,1) - E(1,1)) < tol1(1)) && ...
                    (abs(W(end,2) - E(1,2)) < tol2(1)) && ...
                    (abs(W(end,3) - E(1,3)) < tol3(1))
                p0(:,1) = m;
            else
                p0(:,2) = m;
            end
            mx = (p0(1,1) + p0(1,2))/2;
            my = (p0(2,1) + p0(2,2))/2;
            mz = (p0(3,1) + p0(3,2))/2;
            m = [mx my mz]';
            [T,W] = ode45(@ProgIntegr3D,tspan,m,[],par,options);
            flag1 = 0;
            % If the mid point lies in the basin of attraction of the
            % third stable equilibrium point the routine stops
            if (abs(W(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
                    (abs(W(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
                    (abs(W(end,3) - E(npi+2,3)) < tol3(npi+2))
                flag1 = 1;
                break
            end
            dist =(sqrt((p0(1,1) - p0(1,2))^2 + (p0(2,1) - p0(2,2))^2 + ...
                (p0(3,1) - p0(3,2))^2));
            b = W(end,:); D = [D;b]; M = [M;m'];
        end
        if flag1 == 1
            point1 = [0 0 0]; point2 = [0 0 0]; point3 = [0 0 0];
            return
        end
        % Find in D the index of last trajectory evolving towards the
        % first stable equilibrium
        j1 = 0; j = size(D,1);
        while (j1 == 0)
            if (abs(D(j,1) - E(1,1)) < tol1(1)) && ...
                    (abs(D(j,2) - E(1,2)) < tol2(1)) && ...
                    (abs(D(j,3) - E(1,3)) < tol3(1))
                r = j; j1 = 1;
            else
                j = j - 1;
            end
        end
        % Find in D the index of last trajectory evolving towards the
        % third stable equilibrium
        j1 = 0; j = size(D,1);
        while (j1 == 0)
            if (abs(D(j,1) - E(1,1)) > tol1(1)) || ...
                    (abs(D(j,2) - E(1,2)) > tol2(1)) || ...
                    (abs(D(j,3) - E(1,3)) > tol3(1))
                s = j; j1 = 1;
            else
                j = j - 1;
            end
        end
        % Declare the separatrix point
        point2(1) = (M(r,1) + M(s,1))/2; point2(2) = (M(r,2) + M(s,2))/2;
        point2(3) = (M(r,3) + M(s,3))/2;
        point1(1) = 0; point1(2) = 0; point1(3) = 0;
        point3(1) = 0; point3(2) = 0; point3(3) = 0;
        return
    end
elseif ((abs(Z1(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
        (abs(Z1(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
        (abs(Z1(end,3) - E(npi+2,3)) < tol3(npi+2))) || ...
        ((abs(Z2(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
        (abs(Z2(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
        (abs(Z2(end,3) - E(npi+2,3)) < tol3(npi+2)))
    % Check if the first initial condition evolves towards the second 
    % stable equilibrium point
    if (abs(Z1(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
            (abs(Z1(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
            (abs(Z1(end,3) - E(npi+2,3)) < tol3(npi+2))
        p0 = p1;
    else
        p0(1,1) = p1(1,2); p0(1,2) = p1(1,1); p0(2,1) = p1(2,2); 
        p0(2,2) = p1(2,1); p0(3,1) = p1(3,2); p0(3,2) = p1(3,1);
        p0 = [p0(1,1) p0(1,2); p0(2,1) p0(2,2); p0(3,1) p0(3,2)];
    end
    % Integrate the system with the new initial conditions
    y0 = p0(:,1); [T,Z1] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
    y0 = p0(:,2); [T,Z2] = ode45(@ProgIntegr3D,tspan,y0,[],par,options);
    % Check if the initial conditions evolve towards two different 
    % stable equilibrium points
    if (abs(Z1(end,1) - Z2(end,1)) < tol4b(1))  && ...
            (abs(Z1(end,2) - Z2(end,2)) < tol4b(2)) && ...
            (abs(Z1(end,3) - Z2(end,3)) < tol4b(3))
        point1 = [0 0 0]; point2 = [0 0 0]; point3 = [0 0 0];
        return
    end
    % Necessarly the the second initial condition evolves towards the 
    % third stable equilibrium point
    % Initialize the matrix of trajectories
    D = [Z1(end,:);Z2(end,:)];
    % Initialize the matrix of mid points
    M = [p0(:,1)';p0(:,2)'];
    % Compute the first mid point
    mx = (p0(1,1) + p0(1,2))/2; my = (p0(2,1) + p0(2,2))/2;
    mz = (p0(3,1) + p0(3,2))/2;
    m = [mx my mz]';
    % Integrate the system with the mid point as initial condition
    [T,W] = ode45(@ProgIntegr3D,tspan,m,[],par,options);
    flag1 = 0;
    if nsp == 3
        % If the mid point lies in the basin of attraction of the first
        % stable equilibrium point the routine stops
        if (abs(W(end,1) - E(1,1)) < tol1(1)) && ...
                (abs(W(end,2) - E(1,2)) < tol2(1)) ...
                && (abs(W(end,3) - E(1,3)) < tol3(1))
            point1 = [0 0 0]; point2 = [0 0 0]; point3 = [0 0 0];
            return
        end
        % Update matrices
        D = [D; W(end,:)]; M = [M;m'];
        dist = sqrt((p0(1,1) - p0(1,2))^2+(p0(2,1) - p0(2,2))^2+ ...
            (p0(3,1) - p0(3,2))^2);
        % Bisection routine
        while (dist > tol)
            if (abs(W(end,1) - E(npi+2,1)) < tol1(npi+2)) && ...
                    (abs(W(end,2) - E(npi+2,2)) < tol2(npi+2)) && ...
                    (abs(W(end,3) - E(npi+2,3)) < tol3(npi+2))
                p0(:,1) = m;
            else
                p0(:,2) = m;
            end
            mx = (p0(1,1)+ p0(1,2))/2; my = (p0(2,1)+ p0(2,2))/2;
            mz = (p0(3,1)+ p0(3,2))/2;
            m = [mx my mz]';
            [T,W] = ode45(@ProgIntegr3D,tspan,m,[],par,options);
            if (abs(W(end,1) - E(1,1)) < tol1(1)) && ...
                    (abs(W(end,2) - E(1,2)) < tol2(1)) && ...
                    (abs(W(end,3) - E(1,3)) < tol3(1))
                flag1 = 1;
                break
            end
            dist =(sqrt((p0(1,1) - p0(1,2))^2+(p0(2,1) - p0(2,2))^2+ ...
                (p0(3,1) - p0(3,2))^2));
            b = W(end,:); D = [D;b]; M = [M;m'];
        end
        if flag1 == 1
            point1 = [0 0 0]; point2 = [0 0 0]; point3 = [0 0 0];
            return
        end
        % Find in D the index of last trajectory evolving towards the
        % second stable equilibrium
        j1 = 0; j = size(D,1);
        while (j1 == 0)
            if (abs(D(j,1) - E(npi+2,1)) < tol1(npi+2)) && ...
                    (abs(D(j,2) - E(npi+2,2)) < tol2(npi+2)) && ...
                    (abs(D(j,3) - E(npi+2,3)) < tol3(npi+2))
                r = j; j1 = 1;
            else
                j = j - 1;
            end
        end
        % Find in D the index of last trajectory evolving towards the
        % third stable equilibrium
        j1 = 0; j = size(D,1);
        while (j1 == 0)
            if (abs(D(j,1) - E(npi+2,1)) > tol1(npi+2)) || ...
                    (abs(D(j,2) - E(npi+2,2)) > tol2(npi+2)) || ...
                    (abs(D(j,3) - E(npi+2,3)) > tol3(npi+2))
                s = j; j1 = 1;
            else
                j=j - 1;
            end
        end
        % Declare the separatrix point
        point3(1) = (M(r,1)+M(s,1))/2; point3(2) = (M(r,2)+M(s,2))/2;
        point3(3) = (M(r,3)+M(s,3))/2;
        point2(1) = 0; point2(2) = 0; point2(3) = 0;
        point1(1) = 0; point1(2) = 0; point1(3) = 0;
        return
    else
        point3(1) = 0; point3(2) = 0; point3(3) = 0;
        point2(1) = 0; point2(2) = 0; point2(3) = 0;
        point1(1) = 0; point1(2) = 0; point1(3) = 0;
        return
    end
end