Solving $y' = 2xy^2$ using backward Euler

Question

$$ y' = 2xy^2, \quad y(0) = 1$$

The formula for backward Euler is

$$ y_{n + 1} = y_n + hf(x_{n+1}, y_{n+1}) $$

Where $h$ is the step size.

Plugging in for $f$

$$ y_{n + 1} = y_n + h \times 2 \times x_{n+1} \times y_{n+1}^2 $$

Using the quadratic formula for to solve for $y_{n + 1}$

$$ y_{n + 1} = \frac{1\pm\sqrt{1-4 \times 2 \times h \times (x + h) \times y_n}}{2 \times h \times 2 \times (x + h)} $$

I got the same answer as Wolfram alpha only when I used the negative sign in the quadratic. Why is this?

Here's the Matlab code so you can see for yourself

clc;
clearvars;

h = 0.05;
x = 0;
y = 1;

disp('      x_n     |     y_n     ')
while x <= 2 + h;

   yplus1 = (1 - sqrt(1 - 4*2*h*(x + h)*y))/(2*h*2*(x+h));

   fprintf('%e | %e \n', (x + h), yplus1);

   x = x + h;

   y = yplus1;

end

Wolfram Alpha link

http://www.wolframalpha.com/input/?i=%7By%27%28x%29+%3D+2xy%5E2%2C+y%280%29%3D1%7D+from+0+to+0.5+by+backward+Euler

Answer 1

You are to solve $$ y_{n + 1} = y_n + h·2·x_{n+1}· y_{n+1}^2 $$ where the quadratic term is small and you expect the solution to be close to $y_n$. Multiply with $8hx$ to get $$ 0=8hx·y_n-1+(4hx·y_{n+1}-1)^2 $$ $h$ is supposed to be very small, so $$ \sqrt{1-8hx·y_n}=1-4hx·y_{n+1} \\\iff 4hx·y_{n+1}=1-\sqrt{1-8hx·y_n}=\frac{8hx·y_n}{1+\sqrt{1-8hx·y_n}} $$ and thus finally $$ y_{n+1}=\frac{2·y_n}{1+\sqrt{1-8·h·x_{n+1}·y_n}} $$