I try to understand an old Fortran code that is not well documented. In this code the ODE
$$ \frac{dy}{dx} = -\frac{B(x - y)}{y} $$
is solved numerically as an initial value problem from $x_0=0.99$, $y_0=0.01$ with a step size $h>0$ in negative x-direction. The documentation says that it uses central difference quotient. This would be:
$$ \frac{dy}{dx}=\frac{y(x + h) - y(x - h)}{2h}. $$
The code calculates $y_n=y(x-h)$ iteratively from following equation
$$ y_n = -\left(\frac{B h}{2}\right) + \sqrt{\left(\frac{B h}{2}\right)^2 + B h (2x - y) + y^2} $$
with $y=y(x)$.
I transformed this to:
$$ \begin{split}\Leftrightarrow \left(y_n + \frac{B h}{2}\right)^2 & = \left(\frac{B h}{2}\right)^2 + B h (2x - y) + y^2 \\ \Leftrightarrow y_n^2 + B h y_n - y^2 & = B h (2x - y) \\ \Leftrightarrow y_n^2 + B h y_n - y^2 - B h x & = B h (x - y) \\ \Leftrightarrow \frac{y^2 - y_n^2 + B h(x - y_n)}{y h} & = -\frac{B (x - y)}{y}.\end{split} $$
I would like to know how the left hand side in the last line of these equations relates to the central difference quotient. Is there some approximation used? I also solved this ODE with a forward difference quotient. The results are similar, so I think the formula in the code is correct, but I would like to understand why.