NDSolve giving a wrong solution

The general problem

In using NDSolve to solve first-order IVPs, there are basically two ways to set up the ODE:

y'[x] == f[x, y[x]]     (* explicit form *)
F[x, y[x], y'[x]] == 0  (* implicit form *)

Most numerical solvers require the problem to be specified in explicit form. In Mathematica, there is only one solver that works with the implicit form, IDA, and it is restricted to machine precision. Since it is easy to convert the implicit form to an explicit second-order ODE by differentiating with respect to x, perhaps there has not been much pressure to develop implicit-form solvers.

In Mathematica, you can request that a solution be attempted in either form with the Method option:

Method -> {"EquationSimplification" -> "Solve"}    (* explicit *)
Method -> {"EquationSimplification" -> "Residual"} (* implicit *)

With the "Solve" method, which is the default, NDSolve calls Solve to convert an ODE to explicit form. An equation given in implicit form may have multiple solutions, and if so, NDSolve will integrate each one separately. That is what happens in the OP's example. Further, NDSolve is set up to integrate the separate explict-form ODEs independently and cannot combine them, which is what is required in the OP's case (see @BobHanlon's answer).

Now Solve's issue of genericity plays an important role here. In the OP's case, it returns solutions that are each valid over certain domains and invalid in other nonempty regions, including ones we want to integrate over. Reduce is much more careful and correctly analyzes the OP's system. One can make Solve use Reduce with the option Method -> Reduce, but it still returns two separate solutions, each valid one side of x == 0. Further it returns ConditionalExpression, which NDSolve chokes on (and gives a "non-numerical" NDSolve::ndnum error at the initial condition during the ProcessEquations phase). ConditionalExpression was introduced pretty late, in V8, and perhaps not enough requests to have NDSolve handle it properly have been made to WRI.

OTOH, the "Residual" method solves the ODE implicitly at each step. Since both solutions are valid simultaneously only at x == 0, it will find the right branch once NDSolve takes a step. This computes the correct solution, which the OP mentions. The only drawback is that only one integration method is available and only in machine precision.

It does seem it would be an easy thing in the NDSolve`ProcessEquations stage to check that the original implicit-form ODE is satisfied by the explicit forms at the initial condition. That wouldn't catch the problem in the example at y[0] == 3, at which point both explicit forms satisfy the implicit-form ODE, but it would catch the problem at y[1] == 2. Another issue with the solutions returned by Solve is that the explicit formula for y'[x] needs to switch branches to the other solution returned by Solve when integration crosses x == 0. Switching branches is not something NDSolve is set up to do nor does it seem to me to be an easy programming fix, since each solution is integrated independently. Some ways to do this are given below, but all require the user to prepare the NDSolve call. None are done automatically by NDSolve, which would be desirable.

Finally, what should the user expect? For a long time now in scientific computation, the user has been expected to set up the numerical integration of differential equations. This seems still to be the case in MATLAB and NumPy. I don't know Maple well enough to comment. The general philosophy of Mathematica has been to make everything automatic as much as possible. Mathematica also has tended to use generically true solutions instead of a more rigorous restriction. These are somewhat in conflict here, since the generic methods of Solve are a source of the issue with the NDSolve solutions. On the other hand, to have everything be automatic is not so much a Wolfram goal as a guiding principle. The Q&As on this site show that Automatic does not always get the job done. The user often has to understand the problem, know what solvers are available, prepare the input accordingly, and call the solver with the right options. For an implicit-form IVP, the user should be aware that there can be a problem with solving for y'[x]. They should also be aware that there are standard ways of dealing with implicit-form ODEs:

• using an implicit solver like IDA, called when "Residual" is invoked;
• differentiating to raise the order;
• solving for y'[x] explicitly, the default "Solve" method.

I'll reiterate that I think it's reasonable to expect NDSolve to check that an explicit-form satisfies the original implicit-form ODE at the initial condition. While the user can check the results of NDSolve after the fact, in cases like the IVP y[1] == 2, it would prevent an extraneous integration.

The OP's examples

The explicit solutions for y'[x] we get for the OP's ODE have two branches for x < 0 and two for x > 0. The two solutions result from the (algebraic) rationalization of the ODE, which introduces the possibility of extraneous solutions. In fact, the solution set consists of four connected components, two for the interval x < 0 and two for x > 0. Each solution returned by Solve is valid over one interval but not over both. However, we can transform them into one correct and one incorrect solution by Simplify[..., x > 0], but that's hardly a general technique, I think.

Workaround #1

The OP's discovery:

ode = -x == y'[x] y[x]/(1 + Sqrt[1 + (y'[x])^2]/2);

ListLinePlot[
 NDSolveValue[{ode, y[0] == 3}, y, {x, -7, 7}, 
  Method -> {"EquationSimplification" -> "Residual"}],
 PlotRange -> All
 ]

Workaround #2

Differentiating the ODE raises the order but results in one for which there is a unique explicit form. You have to use the ODE to solve for the initial condition for y'[0].

sol = NDSolve[{D[ode, x], y[0] == 3, y'[0] == 0}, y, {x, -7, 7}]

Workaround #3

Use the correct explicit form, constructed from the correct branches for x <> 0:

ode2 = y'[x] == 
   Piecewise[{
    {(4 x y[x] - Sqrt[3 x^4 + 4 x^2 y[x]^2])/(x^2 - 4 y[x]^2), x < 0}}, 
    (4 x y[x] + Sqrt[3 x^4 + 4 x^2 y[x]^2])/(x^2 - 4 y[x]^2)];

sol = NDSolve[{ode2, y[0] == 3}, y, {x, -7, 7}]

Workaround #4

There are problems with our algebraic notation and its relation to algebraic functions. Applying the assumption x > 0 alters the branch-cut selection when simplifying the solutions returned by Solve so that one of them is correct. In other words, this gives a simpler formula for y'[x] that is equivalent to Workaround #3.

sol = NDSolve[{#, y[0] == 3} /. Rule -> Equal, y, {x, -7, 7}] & /@
  Assuming[x > 0,
   Select[Simplify@Solve[ode, y'[x]], 
    ode /. # /. {y[x] -> 1, x -> 1.`20} &]
   ] // Apply[Join]

Workaround #5

The Solve option Method -> Reduce produces correct solutions in the form of a ConditionalExpression. To get a method that checks and chooses the correct branch of an ODE that implicitly defines y'[x], the user would have to do their own preprocessing. The following is a way in which rhs[] selects the branch that satisfies the original ODE by converting the conditional expressions into a single Piecewise function. The conditions are converted from equations a == b to a comparison Abs[a-b] < 10^-8. I had to add the value at the branch point x == 0 manually.

In other words, this checks y'[x] at each step and selects the correct branch for the step. It will thus automatically switch branches when needed, at x == 0 in the OP's problem. It is worth pointing out that this fixes a problem arising from rationalization of the ODE that introduces extraneous branches. It is possible for an implicit-form ODE to have multiple valid branches. The method below will combine them all (if the solutions have the ConditionalExpression form), which should be considered an error, although it might still accidentally produce a correct solution. For the OP's ODE, it does the right thing.

ClearAll[rhs];
rhs[x_?NumericQ, y_?NumericQ] = Piecewise[
   yp /. Solve[ode /. {y[x] -> y, y'[x] -> yp}, yp, 
       Method -> Reduce] /. ConditionalExpression -> List /. 
    Equal -> (Abs[#1 - #2] < 10^-8 &),
   0 (* y'[0] == 0 *)];

sol = NDSolve[{y'[x] == rhs[x, y[x]], y[0] == 3}, y, {x, -7, 7}]

Here is a very hacky way to fix up the result of the internal Solve result. It is achieved through a sequence of viral UpValues for $tag that rewrites a ConditionalExpression solution into a Piecewise solution like the one above.

opts = Options@Solve;
SetOptions[Solve, Method -> Reduce];
Block[{ConditionalExpression = $tag, $tag},
 $tag /: Rule[v_, $tag[a_, b_]] := $tag[v, a, b];
 $tag /: {$tag[v_, a_, b_]} := $tag[List, v, a, b];
 $tag /: call : {$tag[List, v_, __] ..} := {{v -> 
     Piecewise[
      Unevaluated[call][[All, -2 ;;]] /. $tag -> List /. 
       Equal -> (Abs[#1 - #2] < 1*^-8 &)]}};
 sol = NDSolve[{ode, y[0] == 3}, y, {x, -7, 7}]
 ]
SetOptions[Solve, opts];

How to see what Solve does inside NDSolve

If you want to see what happens internally, you can use Trace. NDSolve uses Solve to solve the ODE for the highest order derivative, if it can, and uses the solution(s) to construct the integral(s). This shows the Solve call and its return value:

Trace[
 NDSolve[
  {ode, y[0] == 3},
  y, {x, -7, 7}],
 _Solve,
 TraceForward -> True,
 TraceInternal -> True
 ]

The reason for that behavior seems to be a big logical bug in NDSolve. During calculation it seems to treat expressions like: y==Sqrt[x] and y^2==x as the same. But, as every user knows here, they're not!

As confirmation, take your particular example: Multiplying by the denominator gives $$-x\left(1-\frac{1}{2} \sqrt{1+(y'(x))^2}\right)=y'(x) y(x).$$ Squaring both sides stupidly and solving for $y'(x)$ creates two branches

NDSolve[{y'[x]==(4 x y[x]+Sqrt[3 x^4 + 4 x^2 y[x]^2])/(x^2 - 4 y[x]^2) , y[0]==3}, y, {x,-6,6}]

and

NDSolve[{y'[x]==(4 x y[x]-Sqrt[3 x^4 + 4 x^2 y[x]^2])/(x^2 - 4 y[x]^2) , y[0]==3}, y, {x,-6,6}]

These are indeed exactly the branches NDSolve does provide, although none is valid.

Even worse, although being fundamental, it does not check the solutions. This would require just an extra line of code in the algorithm as it already uses the tuples $(x_i,y(x_i),y'(x_i)$. Just plug them into the equation and check if its true or false (up to some numerical error).

Edit: NDSolve needs to transform the equation to some kind of standard form, which is controlled by EquationSimplification. There are three possible options for this method: MassMatrix, Residual and Solve which is the default. The latter transforms the equation into a form with no derivatives on one side. The system is then solved with an ordinary differential equation solver. When Residual is chosen all the non-zero terms in the equation are just moved to one side and then solved with a differential algebraic equation solver. This is the reason the result is correct in this case as it doesn't use Solve which is buggy here.

Clear["Global`*"]

sol = DSolve[{-x == y'[x] y[x]/(1 + Sqrt[1 + (y'[x])^2]/2), y[0] == 3}, y, 
   x] // Quiet

(* {{y -> Function[{x}, Sqrt[5 - x^2 + 2 Sqrt[4 - x^2]]]}, 
    {y -> Function[{x}, Sqrt[45 - x^2 - 6 Sqrt[36 - x^2]]]}} *)

FunctionDomain[y[x] /. sol[[1]], x]

(* -2 <= x <= 2 *)

The first solution is valid for -2 <= x <= 2

{-x == y'[x] y[x]/(1 + Sqrt[1 + (y'[x])^2]/2), y[0] == 3} /. sol[[1]] // 
 Simplify[#, -2 <= x <= 2] &

(* {True, True} *)

FunctionDomain[y[x] /. sol[[2]], x]

(* -6 <= x <= 6 *)

The second solution is true for x == 0

{-x == y'[x] y[x]/(1 + Sqrt[1 + (y'[x])^2]/2), y[0] == 3} /. sol[[2]] // 
 FullSimplify[#, -6 <= x <= 6] &

(* {x == 0, True} *)

Plot[Evaluate[y[x] /. sol], {x, -6, 6},
 PlotLegends -> Placed[Automatic, {.75, .2}]]

For the numeric solution, restrict the domain to {- 2, 2}

soln = NDSolve[{-x == y'[x] y[x]/(1 + Sqrt[1 + (y'[x])^2]/2), y[0] == 3}, 
    y, {x, -2, 2}] // Quiet;

The numeric solutions are valid in different portions of the domain

Plot[Evaluate[y[x] /. soln], {x, -2, 2},
 PlotRange -> {0, 3.1},
 PlotLegends -> Placed[Automatic, {.7, .5}]]