function [ET,NB,time_in_state]=repair_facility(N0,T)
t=0;   Nt=N0;  %current number broken between 0 and 4
ET=[0];  NB=[N0];     %initialize event times and number broken vector.
lam=[4 3 2 1 .00000000001];   mu=[.00000000001 2 4 4 4];
while t<T
  TB=-log(rand)/lam(Nt+1);   % time till next breakdown
  TR=-log(rand)/mu(Nt+1);   % time till next repair
  t=t+min(TR,TB);           %update clock
  ET=[ET t];                %update event times
  Nt=Nt+sign(TR-TB);       %number broken after time t.
  NB=[NB Nt];         % update the number broken at time t.
end
d=diff(ET);  n=length(d);  time_in_state=[];
for i=0:4
    time_in_state=[time_in_state sum(d(NB(1:n)==i))];
end
n=length(ET);
x=reshape([(ET(1:(n-1))) ;(ET(2:n))],1,2*(n-1)); y=reshape([(NB(1:(n-1))) ;(NB(1:(n-1)))],1,2*(n-1));
plot(x,y);
xlabel('time') ; ylabel('Number broken')
time_in_state=time_in_state/sum(time_in_state);   % estimates probability of states