function [nstates, Qav, Qre, pi0, av_indx] = CTMC_CR(nm,n,lambdaB,lambdaR,muB,muR)
%CTMC_CR infinitesimal generator matrix and initial probability vector
% of the Component Redundancy (CR) Continuous Time Markov Chain (CTMC)
%   nm       : number of hot standby (n+m>=n+1)
%   n       : number of required rectifiers or batteries to address load
%   lambdaB : battery failure rate
%   lambdaR : rectifier failure rate
%   muB     : battery recovery rate
%   muR     : rectifier recovery rate

aus = nm-n+2;% number of states only considering rectifiers
nstates = (nm-n+2)^2;

eaus = zeros(1,nm-n+2);
eaus(end) = 1;
av_indx = find(kron(ones(1,nm-n+2), eaus)+kron(eaus, ones(1,nm-n+2))<1);

% define rectifiers' and batteries' behavior
Rr_failure = diag(lambdaR*[nm:-1:n],1);
Rb_failure = diag(lambdaB*[nm:-1:n],1);

Rre = kron(Rr_failure, eye(aus, aus)) + kron(eye(aus, aus), Rb_failure);
Rre(av_indx(2:end),1) = muB;

Rav = Rre;
Rav(2:end,1) = muB;

% % Plotting reliability model
% D= digraph(Rre);
% Plot = plot(D);
% highlight(Plot,[5,10,15,20,21,22,23,24,25],'NodeColor','r');
% highlight(Plot,'Edges',[inedges(D,5)',inedges(D,10)',inedges(D,15)',inedges(D,20)',...
%     inedges(D,21)',inedges(D,22)',inedges(D,23)',inedges(D,24)',inedges(D,25)'],'EdgeColor','r');

Qav = Rav - diag(Rav*ones(nstates,1));
Qre = Rre - diag(Rre*ones(nstates,1));

pi0 = zeros(1,nstates);
pi0(1) = 1;

end