(en - m c^2 - e phi ) u == c s ( p - (e/c) A ) v

eqn1 = %

(en + m c^2 - e phi ) v == c s ( p - (e/c) A ) u

eqn2 = %

Eliminate[ {eqn1, eqn2}, v ] ;

Simplify[%]

Solve[%,en] ;

sol = Simplify[%]

Part[Part[Part[sol,2],1],2]

Series[%, {e,0,2} ]

Series[ Expand[ Normal[%]], {p,0,2} ]