clc
Ni=[30886 16548 80194;33607 16117 78730]';
Nm=[47128 26224 139092]';
R =[0.00907 0.0107 0.01233]';
r0=0.029;k=size(Ni);KK=zeros(k);Ai=zeros(k(1),k(2)-1);
for i=1:k(1);
    for j=1:k(2)-1;
    Ai(i,j)=(Ni(i,j)-Ni(i,j)/R(i)*Ni(i,j));
    end
end
A=[ 0             0       Ni(2,1)/Nm(2)       0       Ni(3,1)/Nm(3)       0;
    R(2)/Nm(1)    0       -R(1)/Nm(2)         0          0                 0;
    Ni(1,1)/Nm(1)  0           0              0          0          Ni(3,1)/Nm(3);
    0              0            0        R(3)/Nm(2)        0         -R(2)/Nm(3);
    0            Ni(1,1)/Nm(1)   0       Ni(2,1)/Nm(2)      0               0;
    0           -R(3)/Nm(1)     0           0           R(1)/Nm(3)         0;]
B=[Ai(1,1)-1+Ni(1,1)/Nm(1);     0;
   Ai(2,1)-1+Ni(2,1)/Nm(2);     0;
   Ai(3,1)-1+Ni(3,1)/Nm(3);     0;]
lamd=A^(-1)*B
M=zeros(4,10);M(1,1)=Ni(md1,2);M(2,1)=Ni(2,2);M(3,1)=Ni(3,2);
M(4,1)=Ni(1,2)+Ni(2,2)+Ni(3,2);Nm1=Nm(1);Nm2=Nm(2);Nm3=Nm(3);
for s=1:9
    M(1,s+1)= M(1,s).*(1+R(1).*(1-M(1,s)./Nm1+lamd(3,1).*M(2,s)./Nm2+lamd(5,1).*M(3,s)./Nm3));
    M(2,s+1)= M(2,s).*(1+R(2).*(1-M(2,s)./Nm2+lamd(1,1).*M(1,s)./Nm1+lamd(6,1).*M(3,s)./Nm3));
    M(3,s+1)= M(3,s).*(1+R(3).*(1-M(3,s)./Nm3+lamd(2,1).*M(1,s)./Nm1+lamd(4,1).*M(2,s)./Nm2));
    M(4,s+1)=M(1,s+1)+M(2,s+1)+M(3,s+1);
end
M