
% Parameters of the network
K=3;             % number of jobs
N=2;             % number of nodes
p=[ 0 1; 1 0 ];  % routing probabilities
e=[ 1 1 ];       % visit ratios (by hand)
mu=[ 1 2 ];      % service rates
lambda=[1 1];    % arrival rates at each node
cb=[ 5 0.5 ];    % coefficient of variation of the service times
epsilon=0.01;    % epsilon (termination)
g=0;             % g by Equation 9.17

% STEP 1: Initialization
lambda_l=0;            % lower bound for lambda
lambda_u=min( mu./e ); % upper bound for lambda

% containers for all intermediate results
LAMBDAS=[];   % intermediate values for lambda
LAMBDASLU=[]; % intermediate values for the lower and upper bounds of lambda
KS =[];       % intermediate values for the mean number of jobs
RHOS=[];      % intermediate values for the utilizations
TS=[];        % intermediate values for the mean response time

% STEP 2
while abs( g - K ) > epsilon

  % STEP 2.1: Bisection
  lambda= ( lambda_l + lambda_u ) /2;

  % STEP 2.2

  % The throughputs for each node
  lambdas=e .* lambda;

  % The utilizations
  rho= lambdas./mu;


  cij=[ 1 1 ];
  ca=cij;
  delta_ca=[ ca .+ 2 * epsilon ];
  while min( abs( delta_ca .- ca ) ) > epsilon
    % Phase 1: Merging
    delta_ca=ca;
    ca=cij; % In our case this is quite easy
    % Phase 2: Flow
    cd= rho.^2 .* ( cb .+ 1 ) + ( 1 .- rho ) .* ca + rho .* ( 1 - 2.*rho );
    % Phase 3: Splitting
    cij=cd(end:-1:1);
    %ca=cij
  end

  % Correction parameter (improved)
  a=( cb + ca ) /2;

  % Correction parameter (standard)
  a=( cb + 1 ) /2;

  % System equation
  f=rho .+ ( ( rho.^2 .* a ) ) ./ ( 1 - ( ( K - 1 - a ) / ( K - 1 ) ) .* rho );

  g=sum(f);

  LAMBDASLU=[ LAMBDASLU; lambda_l, lambda_u ];
  KS=[ KS; f ];
  TS=[ TS; f./lambda ];
  RHOS=[ RHOS; rho ];

  % Calculate the new bounds
  if g < K
    lambda_l = lambda;
  else
    lambda_u = lambda;
  end

  LAMBDAS=[LAMBDAS; lambda ];

end
LAMBDAS
LAMBDASLU
KS 
TS 
RHOS
