My first inclination was to stay symbolic all the way through and then compute numerical values at the end.
syms lambda mu sigma positive
syms n k integer
a_n = symprod(lambda/(mu + k*sigma),k,1,n);
a_n = simplify(a_n,10)
P_0 = 1/(1 + symsum(a_n,n,1,inf));
P_0 = simplify(P_0,10)
P_n = P_0*a_n
P_mu = symsum(P_n*(mu/(mu + n*sigma)),n,0,inf)
P_c = (1 - P_0)*(mu/(mu + sigma))
Now evaluate using subs. For example
muval = (1:40).';
Pc_record = subs(P_c, [lambda sigma mu] , {20, 0.1 , muval});
Pmu_record = subs(P_mu, [lambda sigma mu] , {20, 0.1 , muval});
P0_record = subs(P_0, [lambda sigma mu] , {20, 0.1 , muval});
VPA to see the values. Need to increase digits to deal with very large terms in the expressions.
digits(100);
figure
h1 = subplot(211);hold on;plot(muval,real(vpa(Pc_record))),ylabel('real(Pc)')
h2 = subplot(212);hold on;plot(muval,imag(vpa(Pc_record))),ylabel('imag(Pc)');xlabel('mu')
figure
h3 = subplot(211);hold on;plot(muval,real(vpa(Pmu_record))),ylabel('real(Pmu)')
h4 = subplot(212);hold on;plot(muval,imag(vpa(Pmu_record))),ylabel('imag(Pmu)');xlabel('mu')
figure
h5 = subplot(211);hold on;plot(muval,real(vpa(P0_record))),ylabel('real(P0)')
h6 = subplot(212);hold on;plot(muval,imag(vpa(P0_record))),ylabel('imag(P0)'),xlabel('mu')
It looks like for large vaues of mu some error is creeping into the vpa calculations yielding a very, very, small imaginary component even with digits(100). I assume that imaginary part can be safely ignored. Or maybe use more digits? Or maybe the expressions can be further simplified symbolically before subbing in specific values. Sometimes that helps.
Alternatively, setting lambda and sigma ahead of the loop seems to work better, but slower. Couldn't run for all values of mu because it takes too long (takes much longer to run here than on my local machine)
digits(32) % reset to default
clear P0_record Pmu_record Pc_record
syms lambda mu sigma positive
syms n k integer
lambda = sym(20);
sigma = sym(0.1);
Can only use a few values of muval to stay under the Answers runtime limit
muval = [1:8:40 40];
for ii = 1:numel(muval)
mu = sym(muval(ii));
a_n = symprod(lambda/(mu + k*sigma),k,1,n);
a_n = simplify(a_n,10);
P_0 = 1/(1 + symsum(a_n,n,1,inf));
P_0 = simplify(P_0);
P_n = P_0*a_n;
P_mu = symsum(P_n*(mu/(mu + n*sigma)),n,0,inf);
P_c = (1 - P_0)*(mu/(mu + sigma));
Pc_record(ii) = double(P_c);
Pmu_record(ii) = double(P_mu);
P0_record(ii) = double(P_0);
end
plot(h1,muval,real(Pc_record),'-x'),grid
plot(h2,muval,imag(Pc_record),'-x'),grid
plot(h3,muval,real(Pmu_record),'-x'),grid
plot(h4,muval,imag(Pmu_record),'-x'),grid
plot(h5,muval,real(P0_record),'-x'),grid
plot(h6,muval,imag(P0_record),'-x'),grid
