levinson bromberg 2008 eqn 11

import sympy
sympy.init_printing()
b = sympy.Symbol('beta')
g = sympy.Symbol('Gamma')
m_p = sympy.Symbol('m_p')
m_e = sympy.Symbol('m_e')
J = sympy.Symbol('J')
c = sympy.Symbol('c')
ntb = sympy.Symbol("{n'}_b")
npm = sympy.Symbol('n_{\pm}')
u0 = sympy.Symbol('u^0')
ux = sympy.Symbol('u^x')
C1 = sympy.Symbol('C_1')
C2 = sympy.Symbol('C_2')
Tb0x = sympy.Symbol('T_b^{0 x}')
Tpm0x = sympy.Symbol('T_{\pm}^{0 x}')
Tr0x = sympy.Symbol('T_r^{0 x}')
Tbxx = sympy.Symbol('T_b^{x x}')
Tpmxx = sympy.Symbol('T_{\pm}^{x x}')
Trxx = sympy.Symbol('T_r^{x x}')
Ttr00 = sympy.Symbol("{T'}_r^{0 0}")
Ttr0x = sympy.Symbol("{T'}_r^{0 x}")
Ttrxx = sympy.Symbol("{T'}_r^{x x}")
Tr00 = sympy.Symbol('T_r^{0 0}')
xi = sympy.Symbol('xi')
eqn1 = J - m_p*c**2*ntb*ux
eqn2a = (C1 - Tb0x - Tpm0x - Tr0x).subs(Tb0x,ntb*m_p*c**2*u0*ux).subs(Tpm0x,npm*m_e*c**2*u0*ux)
eqn2b = (C2 - Tbxx - Tpmxx - Trxx).subs(Tbxx,ntb*m_p*c**2*ux*ux).subs(Tpmxx,npm*m_e*c**2*ux*ux)
eos = Ttrxx - Ttr00*(1+xi)/3
radiation_stress_energy_lab = sympy.Matrix([[Tr00,Tr0x],[Tr0x,Trxx]])
radiation_stress_energy_rest = sympy.Matrix([[Ttr00,Ttr0x],[Ttr0x,Ttrxx]])
lorentz_boost = sympy.Matrix([[g,g*b],[g*b,g]]).subs(b,-b)
boost_eqn_1 = (radiation_stress_energy_rest-lorentz_boost*radiation_stress_energy_lab*lorentz_boost)[0,0]
boost_eqn_2 = (radiation_stress_energy_rest-lorentz_boost*radiation_stress_energy_lab*lorentz_boost)[0,1]
boost_eqn_3 = (radiation_stress_energy_rest-lorentz_boost*radiation_stress_energy_lab*lorentz_boost)[1,1]
temp = sympy.solve([boost_eqn_2,boost_eqn_1,eos,boost_eqn_3,eqn2b,eqn2a,eqn1],[ntb,Ttr0x,Tr00,Ttr00,Ttrxx,Trxx,Tr0x])[Ttr0x]
temp = temp.subs(u0,g).subs(ux,g*b)
temp.simplify()