%% CODE DESCRIPTION

% This code solves the linear system G tau = -1 to compute the mean exit
% time to reach the stable fixed point starting from the unstable ones for
% the Markov process (macro model), for different sistem size n

% List of parameter
% n: number of nodes
% N: number of possible edges
% c_1: birth rate
% c_2: death rate
% c_3: traidic closure constant
% c_3_tilde: size dependent rate

clear all
clc

% Parameter
n=30:10:100; % vector of different sizes
N=n.*(n-1)*0.5; % vcector of possible existing edges
c_1 = 0.025;
c_2 = 0.25;
c_3 = 0.91;

l=length(n);

%% Computation of fix points of the associated cubic equation
rad = roots([-c_3 c_3 -(c_1+c_2) c_1]);
p1 = rad(3); % lower stable level
p2 = rad(2); % unstable level
p3 = rad(1); % higher stable level


%% Initialization
vals1 = zeros(1,l); % solution for N_1-->N_2
vals2 = zeros(1,l); % solution for N_3-->N_2



%% Cicle on dimension n
for k=1:l
    
    c_3_tilde = c_3/(N(k)^2);
    N_1 = floor(p1*N(k)); % low level
    N_2 = floor(p2*N(k)); %unstable level
    N_3 = floor(p3*N(k)); % high level
    
    %% Construction of generator
    dim = N(k)+1;
    A = zeros(dim,dim);
    for i=2:dim-1
        j=i-1;
        u = c_1*(N(k)-j)+c_3_tilde*j*(j-1)*(N(k)-j);
        w = c_2*j;
        A(i,i-1) = w;
        A(i,i) = -w-u;
        A(i,i+1) = u;
    end

    u0 = c_2*N(k);
    A(1,1) = -u0;
    A(1,2) = u0;

    j = dim-1;
    u = c_1*(N(k)-j)+c_3_tilde*j*(j-1)*(N(k)-j);
    w = c_2*j;
    A(dim,dim-1) = w;
    A(dim,dim) = -w-u;

    
    %% Construction of the matrix of linear system
    el = N_2+1;
    Q = [A(1:el-1,1:el-1),A(1:el-1,el+1:end);...
        A(el+1:end,1:el-1),A(el+1:end,el+1:end)]; 
   
  
    b=-ones(dim-1,1);
    
    %% solution of linear systrem
    res = Q\b;
     
    vals1(k) = res(N_1);
    vals2(k) = res(N_3);
    
end

%% Plot in loglog scale


loglog(n,vals1,'black','Linestyle','--', 'Linewidth',3) % N_1-->N_2
hold on
loglog(n,vals2,'black','Linewidth',3) % N_3-->N_2

% Plot the ratio in loglog scale
figure(2)
loglog(n,vals1./vals2,'black','Linewidth', 3)
 