Solve a ODE with an implicit method I'm trying to solve this (autonomous) ODE, with an implicit method.
\begin{align}
y'(x) & = \sqrt{\frac{2-y(x)}{y(x)}}, \quad y(0)=0. \\[10pt]
x & >0
\end{align}
The book gives this hint: "$y'(x)>0$, so $y(x)>0$, for $x>0$".
I would use "backward Euler's " method, $y_{n+1}=y_{n} + kf(y_{n+1})$, where $k$ is the discretization step.
The problem is that when I use Newton's method, I divide for $y(0)=0$, so I don't know how to move. Just by watching the hint, I would take (using MatLab sintax) $y(1)=0, y(2)=\epsilon$, and then I'll use the standard routine of Euler's method. 
Is there another way to solve it?

EDIT: I used the answer of @gimusi, and get this plot: it seems to be okay now 
 A: You should start from $y_0=y(0)=0$ and then evaluate:
$$y_1=y_0+kf(y_1)$$
assuming a trial value $y_1>0$  and then iterating.
with $k=0.1$ it should converge for $y_1\approx 0.29$.
Here are the first values I obtain:

A: Given that you are using Matlab, you can readily solve the equation
$$y_{n+1}=y_n+k\sqrt{\frac{2-y_{n+1}}{y_{n+1}}}$$
because it is equivalent to the cubic
$$y_{n+1}^3-2y_ny_{n+1}^2+\left(y_n^2+k^2\right)y_{n+1}-2k^2=0$$
so you can use the roots function to solve, being careful to chose the smallest real root greater than $y_n$. Also there is an analytical solution for comparison:
$$x=2\sin^{-1}\left(\sqrt{\frac{y}{2}}\right)-\sqrt{2y-y^2}$$
Solve with Matlab
% backward.m -- Solves differential equation via backward Euler method

% Step size
k = 0.01;
x = [0:k:1.6];
y = zeros(size(x));
% Solve with backward Euler
for i = 2:length(y),
    P = [1 -2*y(i-1) (y(i-1)^2+k^2) -2*k^2];
    r = roots(P);
    y(i) = 2;
    for j = 1:length(r),
        if isreal(r(j)) && r(j) > y(i-1)
            y(i) = min(r(j),y(i));
        end
    end
end
% Analytical solution
v = linspace(0,y(end),length(y));
u = 2*asin(sqrt(v/2))-sqrt(v).*sqrt(2-v);
% Plot
plot(u,v,'k-',x,y,'r-');
xlabel('x');
ylabel('y');
title(['Solutions to Differential Equation, k = ' num2str(k)]);
legend('Exact','Backward Euler','Location','northwest');

And plot:

A: $$y'(x)  = \sqrt{\frac{2-y(x)}{y(x)}}, \quad y(0)=0\qquad x >0$$
Obviously the difficulty in numerical calculus comes from the initial point where $y'(x)$ is infinite. So, we have to start from a point $(x_1,y_1)$ close to $x=0$ but not at exactly $x=0$.  An accurate value for $(x_1,y_1)$ can be obtained from the asymptotic expansion of $y(x)$ close to the origin.
Analytically solving the ODE will make this clear :
$$dx=\sqrt{\frac{y}{2-y}}\:dy \quad\to\quad x=\int\sqrt{\frac{y}{2-y}}\:dy +\text{constant}$$
Taking account of the condition $y(0)=0$ :
$$x=\int_0^y\sqrt{\frac{t}{2-t}}\:dt$$
$$x(y)=-\sqrt{2y-y^2}+\cos^{-1}(1-y)$$
The series expansion close to $(0,0)$ is :
$$x(y)= \frac{\sqrt{2}}{3}y^{3/2}+O(y^{5/2})$$
The series expansion for the inverse function is :
$$y(x)=\left(\frac92\right)^{1/3}x^{2/3}+O(x^{5/3})$$
So, a good starting point will be : $\quad y_1=\left(\frac92\right)^{1/3}x_1^{2/3}\quad$ with arbitrary small $x_1$.
For example, $x_1$ can be of the order of magnitude of a few times the discretization step. This arbitrary choice of $\left(x_1\:,\: y_1=\left(\frac92\right)^{1/3}x_1^{2/3} \right)$ will not affect the final result of the computation of $y(x)$.
