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

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);

M