Let me provide some steps of the derivation here to see where the problem had been.
$$\mathcal{E} \Psi_+ = \frac{\sigma \mathbf{P}}{2m} \left(1- \frac{\mathcal{E}}{2mc^2} \right) \sigma \mathbf{P} \Psi_+$$
Let's change the notation a little:
$$\mathcal{E}=E_c-U$$
$$E_c \Psi_+ = \frac{(\sigma \mathbf{P})^2}{2m} \Psi_+ +U \Psi_+ -\frac{E_c}{4m^2 c^2} (\sigma \mathbf{P})^2 \Psi_+ + \frac{1}{4m^2 c^2} ( \sigma \mathbf{P} U \sigma \mathbf{P} )\Psi_+ $$
$$E_c \Psi_+ = \left(1- \frac{(\sigma \mathbf{P})^2}{4m^2 c^2}\right) \left( \frac{(\sigma \mathbf{P})^2}{2m} +U + \frac{1}{4m^2 c^2} ( \sigma \mathbf{P} U \sigma \mathbf{P} ) \right)\Psi_+ $$
$$E_c \Psi_+ = \left( \frac{(\sigma \mathbf{P})^2}{2m} +U - \frac{(\sigma \mathbf{P})^4}{8m^3 c^2}- \frac{(\sigma \mathbf{P})^2 U}{4m^2 c^2} + \frac{1}{4m^2 c^2} ( \sigma \mathbf{P} U \sigma \mathbf{P} ) \right)\Psi_+ $$
It's easy to identify several correct terms from (1) immediately.
One important note: <PERSON> et al derive the expression (1) sans the <PERSON> term for the case when no vector potential is present. Which is probably why I initially got all the extra terms. 818