% maximum number of jobs in the network
K  =50;
e=[1, 1];

% load dependent service rates for terminals and FES node
mu =[
	1/15.0 .* ones(1,K); 
	1.1161, 1.7858, 2.2323, 2.5512 .* ones(1,K-3)];

% initializing (intermediate) result containers
pi1=zeros( K+1, K+1 ); % conditional probabilities for terminal node
pi2=zeros( K+1, K+1 ); % conditional probabilities for FES node
pi1(1,1)=1;            % STEP 1: initialization (terminals) 
pi2(1,1)=1;            % STEP 1: initialization (FES node)
meanT=[ ];
lambda=[];

% Step 2: iteration
for k = 1:K
  % Step 2.1: compute mean response time
  meanT=[
     meanT;
     [ 1 ./ mu(1,1), sum( ( (1:k) ./ mu(2,1:k) ) .* pi2(1:k,k)')]];
  
  % Step 2.2: compute the throughput
  lambda=[
    lambda;
    k ./ sum( e .* meanT( end, : ) ) ];

  % Step 2.3: compute conditional probabilities
  sum1 = 0;
  sum2 = 0;
  for l = k:-1:1
    pi1( l+1, k+1 ) = ( lambda(end) / mu(1,l) ) * pi1( l, k ) * e(1);
    sum1 = pi1( l+1, k+1 ) + sum1;
    pi2( l+1, k+1 ) = ( lambda(end) / mu(2,l) ) * pi2( l, k ) * e(2);
    sum2 = pi2( l+1, k+1 ) + sum2;
  end
  pi1( 1, k+1 ) = 1 - sum1;
  pi2( 1, k+1 ) = 1 - sum2;
end

fprintf('Load-dependent MVA on the reduced network.');
[(10:10:K)', meanT(10:10:K,2)]
fprintf('   = k        = meanT \n');

