%Two devices
P=[];v=[];w=[];W=[];A=[];
n1=1/20;%density of fail flow of device 1  
n2=2/20;%density of fail flow of device 2 
m1=2/20;%density of repare completion flow of device 1 
m2=6/20;%density of repare completion flow of device 2
D1=5;D2=10;%profit
%The matrix 
P=[1-n1-n2     m1        m2             0
    n1       1-m1-n2     0             m2
    n2         0     1-m2-n1           m1
    0          n2        n1          1-m1-m2];
[v,w]=eig(P);W=diag(w)';
for i=1:4;if abs(W(i)-1)<10^(-5);f=i;end;end;
Sa=v(:,f)./sum(v(:,f));%asymptotic states
D=(D1+D2)*Sa(1)+D2*Sa(2)+D1*Sa(3)%total net profit
figure(1);bar(Sa);xlabel('state number');ylabel('probability')%graph of asymptotic states 