restart;
read `Ifproc.map`:
read `Fbe.map`:

        # Basic equation for the chem. potential mud 
        # for y >= 1; z = -mud/y
Eq := series(Fbe(z, 3/2), z, 4) = y^(-3/2)*Zeta(3/2);
        # we choose the substitution  
        # y^(-3/2) -> 1 - x  and z-> w^2
        #  and invert the series
solw := solve(series(Fbe(w^2, 3/2), w, 7) 
        = (1 - x)*Zeta(3/2), w):
        # then
z := series(solw^2, x, 4);  # because  z = w^2 (z = -mu/theta);
z:= convert(z, polynom): # for next evaluation
          # graf of mud=-z*y on y
plot(If(y <= 1, 0, -y*subs(x = 1 - y^(-3/2), z)), y = 0..3);
         # calculation of dimensionless internal energy 
         # per particle Ed, expansion of Fbe(z,5/2) in y
Fy := subs(x = 1 - y^(-3/2), 
           convert(series(Fbe(z, 5/2), x, 4), polynom));
         # using the function If we can express energy for all y
Ed := 3/2/Zeta(3/2)*y^(5/2)*If(y <= 1, Zeta(5/2), Fy):
         #heat capacity is (in units of k)
C := diff(Ed, y):
         #and than we make graph of heat capacity
plot(C, y = 0..3);