# DSolve seems to hang on some d'Alembert first order ode's. What is the cause and is there a workaround?

I noticed that DSolve (version 13.1 on windows 10) seems to hang (waited for too long) on what should be a straight forward solution method for d'Alembert first order odes. Below is just one example.

The question is, what is the cause, where does it happen, and is there any smart workaround? (just knowing the cause of hang/long time will be good enough. Finding work around is extra credit).

## Code

ClearAll[x, y];
ode = y[x] - x == (y'[x])^2*(1 - 2/3*y'[x]);
DSolve[ode, y[x], x]  (*this seems to hang or take too long. Could not wait *)


To help find out where the difficulty could be, I solved this by hand below. The solutions obtained all verify the ode. Also I see no step in the solution which could have cause this hang or very long step as all the steps are direct and involve just solving separable ode that results in the middle of the steps.

## Verification of hand solution

sol1=y->Function[{x},x+1/3]
sol2=y->Function[{x},- 2*(C[1]-x)^(3/2)/3+C[1]]
sol3=y->Function[{x}, 2*(C[1]-x)^(3/2)/3+C[1]]
ode/.sol1
ode/.sol2//Simplify
ode/.sol3//Simplify


## Hand solution

$$\begin{gather*} \boxed{y-x -\left(y^{\prime}\right)^{2} \left(1-\frac{2 y^{\prime}}{3}\right)=0} \end{gather*}$$ This is d'Alembert ODE. Let $$p=y^{\prime}$$ Then the ode becomes \begin{align*} y = x f(p)+g(p) \tag{1} \end{align*} Where \begin{align*} f(p) &= 1 \end{align*} And \begin{align*} g(p) &= p^{2}-\frac{2}{3} p^{3} \end{align*} Eq (1) now becomes \begin{align*} y = x+ p^{2}-\frac{2}{3} p^{3} \tag{1A} \end{align*} Taking derivative of (1A) w.r.t. $$x$$ and remembering that $$p \left(x \right)$$ is a function of $$x$$ gives \begin{align*} p &= 1+\left(-2 p^{2}+2 p\right) \frac{ \mathop{\mathrm{d}p}}{\mathop{\mathrm{d}x}}\tag{2} \end{align*} The singular solution is when $$\frac{ \mathop{\mathrm{d}p}}{\mathop{\mathrm{d}x}}=0$$. This gives From above \begin{align*} p &= 1 \end{align*} Substituting the above in (1A) gives the singular solution as \begin{align*} y&=x +\frac{1}{3} \end{align*} The general solution is when $$\frac{ \mathop{\mathrm{d}p}}{\mathop{\mathrm{d}x}}\neq 0$$. From (2) this results in \begin{align*} \frac{ \mathop{\mathrm{d}p}}{\mathop{\mathrm{d}x}} &=\frac{p -1}{-2 p^{2}+2 p} \end{align*}
Inverting the above gives \begin{align*} \frac{ \mathop{\mathrm{d}x}}{\mathop{\mathrm{d}p}} &=\frac{-2 p^{2}+2 p}{p -1} \end{align*}
Which simplifies to $$\frac{dx}{dp} = -2 p$$ $$x \left(p \right)$$ is now the dependent variable and $$p$$ as the independent variable. This ODE is now solved for $$x$$. This is separable so easy to solve. Its solution is \begin{align*} x= -p^{2}+c_{1} \tag{3} \end{align*} Hence the $$p$$ solutions are \begin{align*} p_{1} &= \sqrt{c_{1}-x}\\ p_{2} &= -\sqrt{c_{1}-x} \end{align*} Substituting each one of the above solution for $$p$$ in (1A) gives the general solution. Substituting $$p_{1}$$ in (1A) gives the first solution $$y=-\frac{2 \left(c_{1}-x \right)^{\frac{3}{2}}}{3}+c_{1}$$ Substituting $$p_{2}$$ in (1A) gives the second solution $$y=\frac{2 \left(c_{1}-x \right)^{\frac{3}{2}}}{3}+c_{1}$$

ClearAll[withTimedIntegrate];
SetAttributes[withTimedIntegrate, HoldFirst];
withTimedIntegrate[code_, tc_] :=
Module[{$$in}, InternalInheritedBlock[{Integrate}, Unprotect[Integrate]; i : Integrate[___] /; ! TrueQ[$$in] :=
Block[{\$in = True},
TimeConstrained[i, tc, Inactivate[i, Integrate]]];
Protect[Integrate];
code]];

ode = y[x] - x == (y'[x])^2*(1 - 2/3*y'[x]);
withTimedIntegrate[DSolve[ode, y, x], 15]
(* output: 0.5MB, three unevaluated Solve[] of equation with
unevaluated integral  *)


I can't decipher what's going on. I looked into the internal working a bit, but it's too complicated. It's not clear whether d'Alembert method of solution is implemented. DSolve seems to get stuck on a very hard integral.

Differentiating the ODE makes it much easier for DSolve; but one has to solve for the extra integration parameter:

ode = y[x] - x == (y'[x])^2*(1 - 2/3*y'[x]);
sol = # /. First@Solve[ode /. #, Cases[#, _C, Infinity]] & /@
DSolve[D[ode, x], y, x]


Solve::svars: Equations may not give solutions for all "solve" variables.

(*
{{y -> Function[{x}, x + 1/3]},
{y -> Function[{x}, -(2/3) (-x + 2 C[1])^(3/2) + 2 C[1]]},
{y -> Function[{x}, 2/3 (-x + 2 C[1])^(3/2) + 2 C[1]]}}
*)

ode /. sol // Simplify

(*  {True, True, True}  *)
`