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);