How can I fix the error?

I want to compute the steady-state probabilities and the service levels of the call center system. The used equations are mentioned in the attached file. I need to fix the code to achieve the same results in the table mentioned in the file.

[Walter Roberson]

Sorry, the ability to execute MATLAB code from an image of the code is not publicly available. I have not even been asked to beta-test the feature to read code from images.

[Noura]

% C1 Program
%The transition of states occurs due to call arrival
%N1=11,N2=16,N3=9,Lambda1=5,Lambda2=4
%Mu1=0.5,Mu2=0.3,Mu3=0.2
N2=input('Enter the number of servers in group 2 : ');
Lambda2=input('Enter the arrival rate of type 2 : ');
Mu2=input('Enter the service rate of group 2 : ');
syms  k2 Rho2
Rho2=((Lambda2)/(Mu2));
L2=symsum(((Rho2)^(k2))/sym(factorial(k2)),k2,0,N2);
P22=((Rho2)^(N2-1))/((factorial(N2-1))*L2); %P(n2=N2-1)
q12=(Lambda2*P22); %q1-2=q3-5=q6-8=q9-11
N1=input('Enter the number of servers in group 1 : ');
Lambda1=input('Enter the arrival rate of type 1 : ');
Mu1=input('Enter the service rate of group 1 : ');
syms  k1 Rho1
Rho1=((Lambda1)/(Mu1));
L1=symsum(Rho1^k1/sym(factorial(k1)),k1,0,N1);
P11=((Rho1)^(N1-1))/((factorial(N1-1))*L1); %P(n1=N1-1)
q13=(Lambda1*P11); %q1-3=q2-5=q4-8=q7-10
q47=Lambda2; %q4-7=q8-10=q11-12
q69=Lambda1; %q6-9=q8-11=q10-12
N3=input('Enter the number of servers in group 3 : ');
Mu3=input('Enter the service rate of group 3 : ');
syms  k4 Rho4
Rho4=((Lambda2)/(Mu3));
L4=symsum(Rho4^k4/sym(factorial(k4)),k4,0,N3);
P44=((Rho4)^(N3-1))/((factorial(N3-1))*L4); %P1(n3=N3-1)
q24=(Lambda2*P44); %q2-4
syms  k3 Rho3
Rho3=((Lambda1)/(Mu3));
L3=symsum(((Rho3)^(k3))/(sym(factorial(k3))),k3,0,N3);
P33=((Rho3)^(N3-1))/((factorial(N3-1))*L3); %P2(n3=N3-1)
q36=(Lambda1*P33); %q3-6
syms  k5 Rho5
Rho5=(Rho3+Rho4);
L5=symsum(((Rho5)^(k5))/(sym(factorial(k5))),k5,0,N3);
P55=((Rho5)^(N3-1))/((factorial(N3-1))*L5); %P3(n3=N3-1)
Lambda3=(Lambda1+Lambda2);
q58=(Lambda3*P55); %q5-8
%The transition of states occurs due to service fulfilment
q21=(N2*Mu2); %q2-1=q5-3=q8-6=q11-9
q31=(N1*Mu1); %q3-1=q5-2=q8-4=q10-7
q42=(N3*Mu3); %q4-2=q6-3=q8-5
syms  k6 
L6=symsum(((Rho2)^(k6))/(sym(factorial(k6))),k6,0,N2-1);
a1=(N2*((Rho2)^(N2)))/((factorial(N2))*(N2-Rho2));
P02=((L6+a1)^(-1));
a2=(((Rho2)^(N2+1))/(N2*(factorial(N2))));
P66=(a2*P02); %P(n2=N2+1)
q74=((q21+q42)*P66); %q7-4=q10-8
syms  k7 
L7=symsum(((Rho1)^(k7))/(sym(factorial(k7))),k7,0,N1-1);
a3=(N1*((Rho1)^(N1)))/((factorial(N1))*(N1-Rho1));
P01=((L7+a3)^(-1)); 
a4=(((Rho1)^(N1+1))/(N1*(factorial(N1))));
P77=(a4*P01); %P(n1=N1+1)
q96=((q31+q42)*P77); %q9-6=q11-8
q1210=((q31+(0.5*q42))*P77); %q12-10
q1211=((q21+(0.5*q42))*P66); %q12-11
%The equations for the steady-state probabilities
%P1(q1-2+q1-3)=P2q2-1+P3q3-1 (Eqn.1)
b1=(q12+q13);
%P2(q2-1+q2-4+q2-5)=P1q1-2+P4q4-2+P5q5-2 (Eqn.2)
b2=(q21+q24+q13); %q25=q13
%P3(q3-1+q3-5+q3-6)=P1q1-3+P5q5-3+P6q6-3 (Eqn.3)
b3=(q31+q12+q36); %q35=q12
%P4(q4-2+q4-7+q4-8)=P2q2-4+P7q7-4+P8q8-4 (Eqn.4)
b4=(q42+q47+q13); %q48=q13
%P5(q5-2+q5-3+q5-8)=P2q2-5+P3q3-5+P8q8-5 (Eqn.5)
b5=(q31+q21+q58); %q52=q31,q53=q21
%P6(q6-3+q6-8+q6-9)=P3q3-6+P8q8-6+P9q9-6 (Eqn.6)
b6=(q42+q12+q69); %q63=q42,q68=q12
%P7(q7-4+q7-10)=P4q4-7+P10q10-7 (Eqn.7)
b7=(q74+q13); %q710=q13
%P8(q8-4+q8-5+q8-6+q8-10+q8-11)=P4q4-8+P5q5-8+P6q6-8+P10q10-8+P11q11-8 (Eqn.8)
b8=(q31+q42+q21+q47+q69); %q84=q31,q85=q42,q86=q21,q810=q47,q811=q69
%P9(q9-6+q9-11)=P6q6-9+P11q11-9 (Eqn.9)
b9=(q96+q12); %q911=q12
%P10(q10-7+q10-8+q10-12)=P7q7-10+P8q8-10+P12q12-10 (Eqn.10) 
b10=(q31+q74+q69); %q107=q31,q108=q74,q1012=q69
%P11(q11-8+q11-9+q11-12)=P8q8-11+P9q9-11+P12q12-11 (Eqn.11)
b11=(q96+q21+q47); %q118=q96,q119=q21,q1112=q47
%P12(q12-10+q12-11)=P10q10-12+P11q11-12 (Eqn.12)
b12=(q1210+q1211); 
A=[b1 -q21 -q31 0 0 0 0 0 0 0 0 0 0; -q12 b2 0 -q42 -q31 0 0 0 0 0 0 0 0; -q13 0 b3 0 -q21 -q42 0 0 0 0 0 0 0; 0 -q24 0 b4 0 0 -q74 -q31 0 0 0 0 0; 0 -q13 -q12 0 b5 0 0 -q42 0 0 0 0 0; 0 0 -q36 0 0 b6 0 -q21 -q96 0 0 0 0; 0 0 0 -q47 0 0 b7 0 0 -q31 0 0 0; 0 0 0 -q13 -q58 -q12 0 b8 0 -q74 -q96 0 0; 0 0 0 0 0 -q69 0 0 b9 0 -q21 0 0; 0 0 0 0 0 0 -q13 -q47 0 b10 0 -q1210 0; 0 0 0 0 0 0 0 -q69 -q12 0 b11 -q1211 0; 0 0 0 0 0 0 0 0 0 -q69 -q47 b12 0; 1 1 1 1 1 1 1 1 1 1 1 1 1];             
AA=rref(A);
AAA=double(AA);
disp('AAA is : ')
disp(AAA)
AAAA=sum(AAA);
disp('AAAA is : ')
disp(AAAA)
%The service levels
%P1ns=((P9+P11)(sum(i=K1:inf)P(n1=i)))+((P12)(sum(i=K2:inf)P(n1=i))),I=(sum(i=K1:inf)P(n1=i)),II=(sum(i=K2:inf)P(n1=i))
T1=input('Enter the fixed time for which a first-type call service is possible : '); %T1=20
k8=(N1+N3+(T1*(q31+q42))); %K8=K1
k8=round(k8);
k9=(N1+(0.5*N3)+(T1*(q31+(0.5*q42)))); %K9=K2
k9=round(k9);
%I=sum(i=K1:inf)P(n1=i)
%P(n1=i)=(((Rho1)^(i))/((N1!)((N1)^(i-N1)))P0,P0=((((Rho1)^(N1))/(N1!))*((N1)/(N1-Rho1))+(sum(i=0:N1-1)(((Rho1)^(i))/i!)))^(-1)
c2=(((N1)*((Rho1)^(N1)))/((factorial(N1))*(N1-Rho1)));
syms k10 
L8=symsum(((Rho1)^(k10))/sym(factorial(k10)),k10,0,N1-1);
L8=double(L8);
P0=(c2+L8)^(-1);
syms k11
L9=symsum(((Rho1)^(k11))/((N1)^(k11-N1)),k11,k8,Inf);
I=((P0/(factorial(N1)))*L9);
I=double(I);
syms k12
L10=symsum(((Rho1)^(k12))/((N1)^(k12-N1)),k12,k9,Inf);
II=((P0/(factorial(N1)))*L10);
II=double(II);
d1=AAA(9,13);
d2=AAA(11,13);
d3=AAA(12,13);
d4=d1+d2;
P1ns=((d4*I)+(d3*II));
disp('P1ns is : ')
disp(P1ns)
P1sL=1-P1ns;
disp('P1sL is : ')
disp(P1sL)
%P2ns=((P7+P10)(sum(i=K3:inf)P(n2=i)))+((P12)(sum(i=K4:inf)P(n2=i))),III=(sum(i=K3:inf)P(n2=i)),IIII=(sum(i=K4:inf)P(n2=i))
T2=input('Enter the fixed time for which a second-type call service is possible : '); %T2=30
k13=(N2+N3+(T2*(q21+q42))); %K13=K3
k13=round(k13);
k14=(N2+(0.5*N3)+(T2*(q21+(0.5*q42)))); %K14=K4
k14=round(k14);
%III=sum(i=K3:inf)P(n2=i)
%P(n2=i)=(((Rho2)^(i))/((N2!)((N2)^(i-N2)))P0,P0=((((Rho2)^(N2))/(N2!))*((N2)/(N2-Rho2))+(sum(i=0:N2-1)(((Rho2)^(i))/i!)))^(-1)
c3=(((N2)*((Rho2)^(N2)))/((factorial(N2))*(N2-Rho2)));
syms k15 
L11=symsum(((Rho2)^(k15))/sym(factorial(k15)),k15,0,N2-1);
L11=double(L11);
P00=(c3+L11)^(-1); %P00=P0
syms k16
L12=symsum(((Rho2)^(k16))/((N2)^(k16-N2)),k16,k13,Inf);
III=((P00/(factorial(N2)))*L12);
III=double(III);
syms k17
L13=symsum(((Rho2)^(k17))/((N2)^(k17-N2)),k17,k14,Inf);
IIII=((P00/(factorial(N2)))*L13);
IIII=double(IIII);
d5=AAA(7,13);
d6=AAA(10,13);
d7=d5+d6;
P2ns=((d7*III)+(d3*IIII));
disp('P2ns is : ')
disp(P2ns)
P2sL=1-P2ns;
disp('P2sL is : ')
disp(P2sL)

[Noura]

Thank you for your response. I wrote the code in the above comment. I want to fix this code to obtain the following results in the mentioned case as P1sl = 0.8075 and p2sl = 0.9551.