clc; clearvars;
A = 100; % traffic intensity
K = 10:10:200;
lK = length(K);
numSim = 80;
compCostErlang = zeros(numSim,lK);
compCostErlangR = zeros(numSim,lK);
for p = 1:lK % to evaluate k
    for s = 1:numSim
        tic; 
        erlang( A,K(p));
        compCostErlang(s,p) = toc;
        tic;
        erlang_r( A,K(p) );
        compCostErlangR(s,p) = toc;
    end   
end

cC_Erlang = mean(compCostErlang);
cC_ErlangR = mean(compCostErlangR);
p = plot(K',cC_Erlang','b--d',K',cC_ErlangR','g:*');
p(1).MarkerSize = 5;
p(1).LineWidth = 1.5;
p(2).LineWidth = 1.5;
grid on
title({'Computational Cost Analysis Erlang Functions';['Given Traffic Intensity = ',num2str(A)]})
xlabel('Number of channels available in the system')
ylabel({'Computational Cost'; '[seconds]'})