Why does DSolve give 1/0 on this first-order Riccati ODE?

Question

This is from a textbook.

Why does DSolve give 1/0 using Mathematica v13.0.1 on Windows 10? Does it happen using other versions? I only have access to 13.0.1 now.

ClearAll[y, mu, x, lambda, a, b]
ode = y'[x] ==   a*Exp[mu*x]*y[x]^2 + lambda*y[x] - a*b^2*Exp[(mu + 2*lambda)*x]
DSolve[ode, y[x], x]

The solution should be according to Maple 2022:

restart;
ode:=diff(y(x),x)=a*exp(mu*x)*y(x)^2+lambda*y(x)-a*b^2*exp((mu+2*lambda)*x)
sol:=dsolve(ode)

gives

y(x) = -b*(_C1*sinh(1/(lambda+mu)*a*b*exp(x*(lambda+mu)))+cosh(1/(lambda+mu)*a*
b*exp(x*(lambda+mu))))/(_C1*cosh(1/(lambda+mu)*a*b*exp(x*(lambda+mu)))+sinh(1/(
lambda+mu)*a*b*exp(x*(lambda+mu))))*exp(x*(lambda+mu)-mu*x)

In Latex:

$$ y \left(x \right) = -\frac{b \left(c_1 \sinh \left(\frac{a b \,{\mathrm e}^{x \left(\lambda +\mu \right)}}{\lambda +\mu}\right)+\cosh \left(\frac{a b \,{\mathrm e}^{x \left(\lambda +\mu \right)}}{\lambda +\mu}\right)\right) {\mathrm e}^{x \left(\lambda +\mu \right)-\mu x}}{c_1 \cosh \left(\frac{a b \,{\mathrm e}^{x \left(\lambda +\mu \right)}}{\lambda +\mu}\right)+\sinh \left(\frac{a b \,{\mathrm e}^{x \left(\lambda +\mu \right)}}{\lambda +\mu}\right)} $$

Answer 1

I think it's a bug. The point of this answer is to show how the bug arises and how DSolve could have avoided it. The problem is a Riccati equation and not particularly hard to solve by hand, despite all the parameters. After some spelunking I discovered the 1/0 issue had to do with the coefficients of the auxiliary 2nd-order linear equation associated with the Riccati equation ode. If you set DSolve`print = Print, then you get an outline of the analysis of the problem by DSolve (clearly this is meant for internal debugging and not for "showing steps" since some steps are missing). We see early on that DSolve seems to think it's a Riccati equation for it solves the associated 2nd-order auxiliary ODE:

Block[{DSolve`print = Print}, DSolve[ode, y, x]]

...

y''[x] == a^2 b^2 E^(mu x + (2 lambda + mu) x) y[x] 
  + (lambda + mu) y'[x]

...

I traced Power[0, -1], the source of the error, with TraceAbove -> True and discovered the divide-by-zero came from the denominator in dsol (discussed below). The denominator is contructed from the coefficients of the constants of integration in the solution to the auxiliary equation. The buggy calls may be inspected with this trace:

Flatten[
  Quiet@Trace[
    DSolve[ode, y, x],
    _Coefficient,
    TraceInternal -> True]
  ][[-2 ;;]]

Here's the outline of how DSolve approached the problem, with an adjustment to Coefficient:

Block[{q0, q1, q2, f, g, h},
 {q0, q1, q2} = Together[CoefficientList[
    -a b^2 E^((2 lambda + mu) x) + lambda y[x] + a E^(mu x) y[x]^2,
    y[x]], Trig -> True];
 dsol2 = DSolve[
    y''[x] == -q0*q2*y[x] + (q1 + D[q2, x]/q2) y'[x],
    y[x], x];
 sol0 = y[x] /. First[dsol2] // ExpandAll;
 (****)
 {f, g, h} = {
   Coefficient[sol0, Sqrt[C[1]]], 
   Coefficient[sol0 /. C[2] -> Log[C[2]], C[2]],
   q2} /. _C -> 1;
 (****)
 dsol = {y[x] -> -((C[1] D[f, x] + D[g, x])/(h (C[1] f + g)))} // 
    Simplify // DSolve`DSolveToPureFunction
 ]

(* InverseFunction::ifun messages omitted...
{y -> Function[{x},
  (b*E^(lambda*x) (E + E^2 C[1] + 
     E^((2 a*b*E^((lambda + mu) x))/(lambda + mu)) C[1]))/
   (E + E^2 C[1] - 
      E^((2 a*b*E^((lambda + mu) x))/(lambda + mu)) C[1])]}
*)

The solution is dsol:

ode /. dsol // Simplify
(*  True  *)

The bug arises in the line highlighted with (****). What DSolve does is to try to get the coefficients of the solution sol0 to the auxiliary equation with the following:

Coefficient[sol0, C[1]]
Coefficient[sol0, C[2]]

But sol0 has the form

sol0
(*
(E^((a b E^(lambda x + mu x))/(lambda + mu) - C[2]) Sqrt[C[1]])/Sqrt[2] - 
 (E^(-((a b E^(lambda x + mu x))/(lambda + mu)) + C[2]) Sqrt[C[1]])/Sqrt[2]
*)

And both calls to Coefficient[] yield zero, since C[1] and C[2] are hidden inside Sqrt[] and E^... respectively. This leads to (C[1] f + g) being zero in the denominator of dsol, which gives the Power::infy error. Constants of integration in linear ODEs (the auxiliary 2nd-order ODE in this case) often end up as an argument of an exponential function, so I would expect some ability of DSolve to handle the. For instance, the following modification yields another way to obtain a solution:

{f, g, h} = {
  sol0 /. {C[1] -> 1, C[2] -> 0}, 
  sol0 /. {C[1] -> 1, C[2] -> 1},
  q2}

One has to avoid C[1] -> 0, so finding a basis for the second-order solution space takes some care. Similarly, one could use SeriesCoefficient or Series. Again the Sqrt[C[1]] poses a slight problem, since the coefficient one wants is for the 1/2 power, not the linear term.

Answer 2

Another way to get a solution, is to solve the ODE for y[x].

ode = y'[x] ==
   a*Exp[mu*x]*y[x]^2 + lambda*y[x] - a*b^2*Exp[(mu + 2*lambda)*x];

rule1 = Thread[CoefficientList[ode[[2]], y[x]] -> {r, s, t}]

(*     {-a b^2 E^((2 lambda + mu) x) -> r, lambda -> s, a E^(mu x) -> t} *)

red = Reduce[
  y'[x] == (ode[[2]] /. rule1) && r != 0 && s != 0 && t != 0, y[x]]

(*   t != 0 && (y[x] == (-s - Sqrt[
     s^2 - 4 r t + 4 t Derivative[1][y][x]])/(2 t) ||
   y[x] == (-s + Sqrt[s^2 - 4 r t + 4 t Derivative[1][y][x]])/(2 t)) &&
  r s != 0   *)

ode2 = red[[2, 1]] /. Reverse /@ rule1

ysol1[c1_] = y /. Flatten@DSolve[ode2, y, x] /. C[1] -> c1

(*   Function[{x}, -((
  E^(-mu x) (I a b E^((lambda + mu) x)
       Cosh[(a b E^((lambda + mu) x))/(lambda + mu)] +
     a b E^((lambda + mu) x)
       c1 Sinh[(a b E^((lambda + mu) x))/(lambda + mu)]))/(
  a (c1 Cosh[(a b E^((lambda + mu) x))/(lambda + mu)] +
     I Sinh[(a b E^((lambda + mu) x))/(lambda + mu)])))]   *)

This satisfies the ODE. The second part of red yields the same solution.

ode /. y -> ysol1[c1] // FullSimplify

(*   True   *)

If you want a real solution:

ceim = ComplexExpand[Im[ysol1[c1][x]] // FullSimplify,
  TargetFunctions -> {Re, Im}]

Solve[ceim == 0, c1]

(*   {{c1 -> 0}}   *)

ysol1[0][x] // FullSimplify

(*   -b E^(lambda x) Coth[(a b E^((lambda + mu) x))/(lambda + mu)]   *)

A first glance shows how to derive this from the solution MichealE2 got.

yME2 = Function[{x}, (b E^(lambda x) (E^C[2] +
         E^((2 a b E^((lambda + mu) x))/(lambda + mu)) Sqrt[C[1]] +
         E^(2 C[2]) Sqrt[C[1]]))/(E^C[2] -
       E^((2 a b E^((lambda + mu) x))/(lambda + mu)) Sqrt[C[1]] +
       E^(2 C[2]) Sqrt[C[1]])] /. {C[1] -> c1, C[2] -> c2};

Solve[yME2[x] == ysol1[0][x], {c1, c2}]

(*   Solve::ifun: Inverse functions are being used by Solve, so some solutions may

{{c1 -> E^(2 c2)/(-1 + E^(2 c2))^2}}   *)

Answer 3

Fyi, this has been fixed in 13.1. Now it does not give 1/0 but gives the correct solution.

The following is screen shot