%-------------------------------------------------------------------------%
%         This software is licensed by Creative Commons BY-NC-SA:         %
%      http://creativecommons.org/licenses/by-nc-sa/3.0/it/legalcode      %
%-------------------------------------------------------------------------%
%
% The file is an example concerning the determination of the basins of 
% attraction of three equilibrium points
%
%-------------------------------------------------------------------------%
clear all
close all
clc 
ngrid = 4; % Number of equispaced points on the edges of the cube
tol = 0.001; % Tolerance used to stop bisection
ed = 6; % Edge of the cube 
ep = [0.1, 0.09, 0.08]; % Shape parameters
npu = [3,4,4]; % Number or PU subdomains
npi = 2; % Number of equilibria partitioning the phase state (tipically
% origin and coexistence equilibrium point)
% Define the system parameters
p = 1; u = 3; a = 5; b = 4; q = 2; v = 2; c = 3; e = 7; r = 2; w = 1;
f = 7; g = 10;
par = [p,u,a,b,q,v,c,e,r,w,f,g]; % Model parameters
% Define the matrix of equilibria
% 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
E = [u 0 0; 
    0 0 0;...
    -u*(p*r*q-g*w*p*v*e-a*v*r*q+a*v*r*w*e+b*w*q*g*v-b*w*r*q)/...
    (-p*r*q+r*c*u*v*a+f*u*w*b*q-f*u*w*a*v*e-g*w*v*c*u*b+g*w*p*v*e)...
    -v*(p*r*q-p*r*c*u-f*u*w*b*q+r*w*c*u*b-r*w*p*e+f*u*w*p*e)/...
    (-p*r*q+r*c*u*v*a+f*u*w*b*q-f*u*w*a*v*e-g*w*v*c*u*b+g*w*p*v*e)...
    (-c*u*p*g*v+r*c*u*v*a+u*p*f*q-u*v*a*f*q+p*q*g*v-p*r*q)*w/...
    (-p*r*q+r*c*u*v*a+f*u*w*b*q-f*u*w*a*v*e-g*w*v*c*u*b+g*w*p*v*e);
    0 v 0; 
    0 0 w;
    u*q*(-p+a*v)/(-p*q+c*u*v*a) p*v*(-q+c*u)/(-p*q+c*u*v*a) 0; 
    u*r*(-p+b*w)/(-p*r+f*u*w*b) 0 w*p*(-r+f*u)/(-p*r+f*u*w*b);...
    0  v*r*(-q+w*e)/(-q*r+g*v*w*e) w*q*(-r+g*v)/(-q*r+g*v*w*e)];
tmax = 90; % Integration time
k = [7,8,6]; % Used to estimate normals
nsp = 3; % Number of stable equilibrium points
% Reconstruct the attraction basins
[A1,A2,A3] = SeparatrixSurfaces(ngrid,tol,ed,ep,npu,npi,par,E,tmax,k,nsp);