# The following statements compute the sojourn time trivariate 
# cumulative distribution function associated with the 
# five cases for the first three customers in an M/M/1 
# queue that is initially empty and idle.  For each case
# the partial derivatives are taken resulting in the 
# trivariate probability distribution functions.  These
# PDFs are integrated over their respective regions to
# calculate their contribution to the trivariate
# distribution.  The contributions are summed for all
# five cases providing the total of 1 as required.

restart; 
assume(lambda > 0); 
assume(mu > 0); 

#  Case A 

FA := simplify((1 - exp(-(lambda + mu) * t1)) * 
     (1 - exp(-(lambda + mu) * t2)) * 
     (1 - exp(-mu * t3))):

#  Case B 
# 
#  FBtop is the joint cdf of t1, t2, t3 when t2 < t3 
#  FBbot is the joint cdf of t1, t2, t3 when t2 > t3 

FBtop := simplify(int( 
            (1 - exp(-(lambda + mu) * t1)) * 
            (1 - exp(-(lambda + mu) * (t2 - y))) * 
            (1 - exp(-mu * (t3 - y))) * 
            mu * exp(-mu * y), 
            y = 0 .. t2)): 

FBbot := simplify(int( 
            (1 - exp(-(lambda + mu) * t2)) * 
            (1 - exp(-(lambda + mu) * (t1 - y))) * 
            (1 - exp(-mu * (t3 - y))) * 
            mu * exp(-mu * y), 
            y = 0 .. t3)):

#  Case C 
# 
#  FCtop is the joint cdf of t1, t2, t3 when t1 < t2 
#  FCbot is the joint cdf of t1, t2, t3 when t1 > t2 

FCtop := simplify(int( 
            (1 - exp(-(lambda + mu) * (t1 - y))) * 
            (1 - exp(-(lambda + mu) * (t2 - y))) * 
            (1 - exp(-mu * t3)) * 
            (lambda + mu) * exp(-(lambda + mu) * y), 
            y = 0 .. t1)): 

FCbot := simplify(int( 
            (1 - exp(-(lambda + mu) * (t1 - y))) * 
            (1 - exp(-(lambda + mu) * (t2 - y))) * 
            (1 - exp(-mu * t3)) * 
            (lambda + mu) * exp(-(lambda + mu) * y), 
            y = 0 .. t2)):

#  Case D 
# 
#  FD1 is the joint cdf of t1, t2, t3 in Case 1 
#  FD2 is the joint cdf of t1, t2, t3 in Case 2 
#  FD3 is the joint cdf of t1, t2, t3 in Case 3 
#  FD4 is the joint cdf of t1, t2, t3 in Case 4 
#  FD5 is the joint cdf of t1, t2, t3 in Case 5 

integrand := (1 - exp(-(lambda + mu) * (t1 - y))) * 
            (1 - exp(-(lambda + mu) * (t2 - y - z))) * 
            (1 - exp(-mu * (t3 - z))) * 
            (lambda + mu) * exp(-(lambda + mu) * y) * 
            mu * exp(-mu * z): 

FD1 := simplify(int(int( 
              integrand, 
              z = 0 .. (t2 - y)), 
              y = 0 .. t2)): 

FD2 := simplify(int(int( 
              integrand, 
              z = 0 .. (t2 - y)), 
              y = 0 .. t1)): 

FD3 := simplify(int(int( 
              integrand, 
              y = 0 .. (t2 - z)), 
              z = 0 .. t3)):

FD4 := simplify(int(int( 
              integrand, 
              z = 0 .. t3), 
              y = 0 .. (t2 - t3)) + 
      int(int( 
              integrand, 
              z = 0 .. (t2 - y)), 
              y = (t2 - t3) .. t1)):

FD5 := simplify(int(int( 
              integrand, 
              z = 0 .. t3), 
              y = 0 .. t1)): 

#  Case E 
# 
#  FE1 is the joint cdf of t1, t2, t3 in Case 1 
#  FE2 is the joint cdf of t1, t2, t3 in Case 2 
#  FE3 is the joint cdf of t1, t2, t3 in Case 3 
#  FE4 is the joint cdf of t1, t2, t3 in Case 4 
#  FE5 is the joint cdf of t1, t2, t3 in Case 5 

integrand := (1 - exp(-(lambda + mu) * (t1 - y - z))) * 
            (1 - exp(-mu * (t3 - z - w))) * 
            (lambda + mu) * exp(-(lambda + mu) * y) * 
            mu * exp(-mu * z) * 
            mu * exp(-mu * w): 

FE1 := simplify(int(int(int( 
                  integrand, 
                  w = 0 .. (t2 - y - z)), 
                  z = 0 .. (t2 - y)), 
                  y = 0 .. t2)): 

FE2 := simplify(int(int(int( 
                  integrand, 
                  w = 0 .. (t2 - y - z)), 
                  z = 0 .. (t1 - y)), 
                  y = 0 .. t1)): 

FE3 := simplify(int(int(int( 
                  integrand, 
                  y = 0 .. (t2 - w - z)), 
                  w = 0 .. (t3 - z)), 
                  z = 0 .. t3)): 

FE4a := simplify(int(int(int( 
                  integrand, 
                  w = 0 .. (t3 - z)), 
                  z = 0 .. (t1 - y)), 
                  y = 0 .. (t2 - t3)) + 
       int(int(int( 
                  integrand, 
                  w = 0 .. (t2 - y - z)), 
                  z = 0 .. (t1 - y)), 
                  y = (t2 - t3) .. t1)): 

FE4b := simplify(int(int(int( 
                  integrand, 
                  w = 0 .. (t3 - z)), 
                  z = 0 .. (t1 - y)), 
                  y = 0 .. t1)): 


FE5a := simplify(int(int(int( 
                  integrand, 
                  y = 0 .. (t1 - z)), 
                  z = 0 .. (t3 - w)), 
                  w = 0 .. (t2 - t1)) + 
       int(int(int( 
                  integrand, 
                  y = 0 .. (t2 - w - z)), 
                  z = 0 .. (t3 - w)), 
                  w = (t2 - t1) .. t3)): 

FE5b := simplify(int(int(int( 
                  integrand, 
                  y = 0 .. (t1 - z)), 
                  z = 0 .. (t3 - w)), 
                  w = 0 .. t3)): 
> #  Check the PDF in the A case 

fA := simplify(diff(diff(diff(FA, t1), t2), t3)): 
simplify(int(int(int(fA, t1 = 0 .. infinity), t2 = 0 .. infinity), t3 = 0 .. infinity)): 

#  Check the PDF in the B cases 

fBtop := simplify(diff(diff(diff(FBtop, t1), t2), t3)):
fBbot := simplify(diff(diff(diff(FBbot, t1), t2), t3)): 

pBtop := simplify(int(int(int(fBtop, t2 = 0 .. t3), t3 = 0 .. infinity), t1 = 0 .. infinity)): 
pBbot := simplify(int(int(int(fBbot, t3 = 0 .. t2), t2 = 0 .. infinity), t1 = 0 .. infinity)): 
simplify(pBtop + pBbot):

#  Check the PDF in the C cases 

fCtop := simplify(diff(diff(diff(FCtop, t1), t2), t3)): 
fCbot := simplify(diff(diff(diff(FCbot, t1), t2), t3)): 

pCtop := simplify(int(int(int(fCtop, t1 = 0 .. t2), t2 = 0 .. infinity), t3 = 0 .. infinity)): 
pCbot := simplify(int(int(int(fCbot, t2 = 0 .. t1), t1 = 0 .. infinity), t3 = 0 .. infinity)): 
simplify(pCtop + pCbot): 

#  Check the PDF in the D cases 

fD1 := simplify(diff(diff(diff(FD1, t1), t2), t3)): 
fD2 := simplify(diff(diff(diff(FD2, t1), t2), t3)): 
fD3 := simplify(diff(diff(diff(FD3, t1), t2), t3)): 
fD4 := simplify(diff(diff(diff(FD4, t1), t2), t3)): 
fD5 := simplify(diff(diff(diff(FD5, t1), t2), t3)): 

pD1a := simplify(int(int(int(fD1, t2 = 0 .. t1), t1 = 0 .. t3), t3 = 0 .. infinity)): 
pD1b := simplify(int(int(int(fD1, t2 = 0 .. t3), t3 = 0 .. t1), t1 = 0 .. infinity)):
pD1 := simplify(pD1a + pD1b):
pD2 := simplify(int(int(int(fD2, t1 = 0 .. t2), t2 = 0 .. t3), t3 = 0 .. infinity)): 
pD3 := simplify(int(int(int(fD3, t3 = 0 .. t2), t2 = 0 .. t1), t1 = 0 .. infinity)): 
pD4a := simplify(int(int(int(fD4, t2 = t3 .. (t1 + t3)), t1 = 0 .. t3), t3 = 0 .. infinity)): 
pD4b := simplify(int(int(int(fD4, t2 = t1 .. (t1 + t3)), t3 = 0 .. t1), t1 = 0 .. infinity)): 
pD4 := simplify(pD4a + pD4b): 
pD5 := simplify(int(int(int(fD5, t2 = (t1 + t3) .. infinity), t1 = 0 .. infinity), t3 = 0 .. infinity)): 
simplify(pD1 + pD2 + pD3 + pD4 + pD5): 

#  Check the PDF in the E cases 

fE1 := simplify(diff(diff(diff(FE1, t1), t2), t3)): 
fE2 := simplify(diff(diff(diff(FE2, t1), t2), t3)): 
fE3 := simplify(diff(diff(diff(FE3, t1), t2), t3)): 
fE4a:= simplify(diff(diff(diff(FE4a, t1), t2), t3)): 
fE4b:= simplify(diff(diff(diff(FE4b, t1), t2), t3)):
fE5a := simplify(diff(diff(diff(FE5a, t1), t2), t3)): 
fE5b := simplify(diff(diff(diff(FE5b, t1), t2), t3)): 

pE1a := simplify(int(int(int(fE1, t2 = 0 .. t1), t1 = 0 .. t3), t3 = 0 .. infinity)): 
pE1b := simplify(int(int(int(fE1, t2 = 0 .. t3), t3 = 0 .. t1), t1 = 0 .. infinity)): 
pE1 := simplify(pE1a + pE1b): 
pE2 := simplify(int(int(int(fE2, t1 = 0 .. t2), t2 = 0 .. t3), t3 = 0 .. infinity)): 
pE3 := simplify(int(int(int(fE3, t3 = 0 .. t2), t2 = 0 .. t1), t1 = 0 .. infinity)): 
pE4a:= simplify(int(int(int(fE4a, t2 = t3 .. (t1 + t3)), t1 = 0 .. t3), t3 = 0 .. infinity)): 
pE4b:= simplify(int(int(int(fE4b, t2 = (t1 + t3) .. infinity), t1 = 0 .. t3), t3 = 0 .. infinity)): 
pE4 := simplify(pE4a + pE4b):
pE5a := simplify(int(int(int(fE5a, t2 = t1 .. (t3 + t1)), t3 = 0 .. t1), t1 = 0 .. infinity)):
pE5b := simplify(int(int(int(fE5b, t2 = t3 .. (t1 + t3)), t1 = 0 .. t3), t3 = 0 .. infinity)): 
pE5 := simplify(pE5a + pE5b):
simplify(pE1 + pE2 + pE3 + pE4 + pE5):