from scipy.optimize import fsolve
# Variables
R = 8314.3
b = 0.0306 #m**3/kmol
a = 0.548*10**6 #pa m**6/kmol**6
T = 973.1
P = 25*10**6 #Pa
# Calculations
Vi = R*T/P
#x = poly(0,'x')
def f(x):
return P*x**2 *(x-b) +a*(x-b) - R*T*(x**2)
#vec = roots(P*x**2 *(x-b) +a*(x-b) - R*T*(x**2))
#vec = roots([P,(P*b+R*T),a,-a*b])
vec = fsolve(f,0,full_output=1)
volume = vec[0]
dH = 8.0906*10**6 -P*volume +0.548*10**6 /volume
# Results
print "Change in enthalpy = %.2e J/kmol"%(dH)
# Note : Answer is different because fsolve function calculates differently. Please check manually.