import matplotlib.pyplot as grf #Parameters R = 8.31446 mu = 0.02897#mass of mole gamma = 1.4 #Initial Conditions T0 = 273.#chemistry normal conditions P0 = 101325.#760mm V0 = R*T0/P0 ro0 = mu/V0 const = P0*V0**gamma dV = 0.0001 V = V0 P = P0 ro = ro0 dP = const/(V-dV)**gamma - P dro = mu/(V-dV) - ro c = (dP/dro)**0.5 x = V**(1./3.) t = x/c x_ = [x] t_ = [t] t1_ = [t] while V>0.0005: V = V-dV P = R*T0/V constT0 = P*V**gamma P1 = const/V**gamma ro = mu/V dP = constT0/(V-dV)**gamma - P dP1 = const/(V-dV)**gamma - P1 dro = mu/(V-dV) - ro c = (dP/dro)**0.5 c1 = (dP1/dro)**0.5 x = V**(1./3.) t = x/c t1 = x/c1 x_.append(x) t_.append(t) t1_.append(t1) grf.plot(x_,t_,'green') grf.grid() #grf.twinx() grf.plot(x_,t1_,'red') grf.show()