1. Your calculations are correct, except that there is no $m_ec$ term on both sides of momentum conservation (supposedly the electron is at rest both before and after the collision, so its momentum is zero). If it’s not obvious to you that the electron being at rest afterwards implies the smallest possible energy for the photon, give the electron an arbitrary final velocity and prove that it can only increase the energy of the photon (note that you don’t necessarily need to compute an exact answer for this). Other than that, yes, solve a quadratic equation for $E_\gamma$ and you get the exact answer, at least as long as you neglect relativistic corrections.
2. Note, though, that the mass of the nucleus is large compared to the other mass scale in the problem, $E_b/c^2$ ($Mc^2$ is several GeV or more—$m_pc^2$ being little less than a GeV for a single proton—while $E_b$ or $E_\gamma$ are a couple keV at most). We expect or at least feel that the second term in your last equation should be small compared to the first one,
$$
E_\gamma^2/2Mc^2 \ll E_b\,.
$$
OK, so now we have an equation for $E_\gamma$ with the general structure
$$
E_\gamma = f(E_\gamma)\equiv\text{(constant)} +\text{(small term that depends on $E_\gamma$)}\,.
$$
An equation of this form is called a fixed point equation because it says that the sequence $E$, $f(E)$, $f(f(E))$, etc. stays fixed when you pick $E = E_\gamma$.
There is a simple way to arrive at an approximate solution for such an equation: pick an initial approximation $E^{(0)}_\gamma$ and set $E^{(1)}_\gamma = f(E^{(0)}_\gamma)$, $E^{(2)}_\gamma = f(E^{(1)}_\gamma)$, etc. If the $E^{(n)}_\gamma$s converge somewhere, they’ll converge to a solution (because the difference between the LHS and the RHS will become smaller and smaller as $n$ increases); and they will converge if you’re lucky and picked a good initial $E^{(0)}_\gamma$.
In our case, the choice of initial approximation is obvious: set $E^{(0)}_\gamma =\text{(constant)} = E_b$. Then
$$
E^{(1)}_\gamma = E^\mathstrut_b + E_b^2/2Mc^2
$$
is the approximate answer you wanted to obtain, and
$$
E^{(2)}_\gamma - E^{(1)}_\gamma = E_b^3/2(Mc^2)^2 + E_b^4/8(Mc^2)^3\approx E_b^3/2(Mc^2)^2
$$
is a good guess as to the error of that answer.
3. Now all of that was unnecessary in this particular case, of course; we could have just as well solved the quadratic equation and approximated the solution (ignoring the root of order $Mc^2$ that makes no sense with the nonrelativistic equations you used):
\begin{align}
E_\gamma &{}= Mc^2 -\sqrt{(Mc^2)^2 - 2Mc^2E_b}\\
&{}= Mc^2\,\left(1 - [1 - 2E_b/Mc^2]^{1/2}\right)\\
&{}\stackrel{\text{(B)}}{=} Mc^2\,\Bigl(1 -\Bigl[1\\
&\qquad{}+ (1/1!)(1/2)(-2E_b/Mc^2) \\
&\qquad{}+ (1/2!)(1/2)(1/2-1)(-2E_b/Mc^2)^2 \\
&\qquad{}+ (1/3!)(1/2)(1/2-1)(1/2-2)(-2E_b/Mc^2)^3 +\cdots\Bigr]\Bigr)\\
&{}= Mc^2 - Mc^2 + E^\mathstrut_b + E_b^2/2Mc^2 + E_b^3/2(Mc^2)^2 +\cdots\\
&{}= E^\mathstrut_b + E_b^2/2Mc^2 + E_b^3/2(Mc^2)^2 +\cdots
\end{align}
Here, equality (B) uses the binomial formula for exponent $1/2$. Yes, it works for fractional exponents as well, except that you obtain a (potentially) infinite sum instead of a finite one.
4. You see that the result is the same (yay! math works!), but I think successive approximation is clearer. More importantly, it does not require you to find the exact solution, so it may be used when you can’t (which is most of the time). Also, finding an exact solution only to approximate it is a little bit silly (if sometimes necessary): the very point of approximate solutions is that they should be easier to find.
