DSolve gives the wrong solution for an inhomogeneous differential equation

I have the following inhomogeneous differential equation (sorry, I cannot come up with a simpler example)

Hom = -(-1 +
z) z^4 (8 (-5 + 6 z) f0[z] + (20 - 70 z + 52 z^2) f0'[
z] + (-1 + z) z (2 (-5 + 7 z) f0''[z] + (-1 + z) z f0'''[z]))
InHom = -((24 z (24 + z (-48 + z (23 + 2 z))))/(-1 + z)) + (
48 Log[1 - z] (z (18 + z (-24 + 5 z)) + 3 (-2 + z) (-1 + z) Log[1 - z]))/z
myeq = Hom + InHom


DSolve gives me a solution but it doesn't work, meaning that

DSolve[myeq == 0, f0[z], z][[1, 1, 2]];
myeq /. f0 -> ((% /. z -> #) &) // FullSimplify
% === 0


gives False. Any idea why this is happening?

• Could be a bug... but the issue unfortunately seems to depend sensitively on the form of the rather complicated equation you are solving. Commented Nov 26, 2019 at 18:15
• @Bill, the range of $z$ i'm interested is $0<z<1$, I tried to use Assuming but it doesn't change. The solution I'm interested in should be $0$ at $z=0$ and for small $z$ behave linearly in $z$. Commented Nov 26, 2019 at 19:57
• Yeah, it says that the solution diverges as $z^{-5}$ but I don't trust that. For example, I can easily solve the differential equation as a power series. For example Nmax = 15; powser = ((1/z^3 myeq) /. f0 -> ((Sum[\[Alpha][n] #^n, {n, 1, Nmax}]) &)); li = CoefficientList[Series[powser, {z, 0, Nmax}], z] // Normal; lis = Solve[li == 0 li]; fsol = (Sum[\[Alpha][n] z^n, {n, 1, Nmax}]) /. lis[[1]]; Series[myeq /. f0 -> (fsol /. z -> # &), {z, 0, Nmax}] gives the first Nmax terms of the solution I'm interested in. It does not diverge at $z=0$ Commented Nov 26, 2019 at 21:00
• Your checking is wrong. Guess what aaaa === 0 will evaluate to? Commented Nov 27, 2019 at 4:06
• In fact, a simpler example is Hom + InHom[[2]]. Commented Nov 27, 2019 at 5:27

With Hom, InHom, and myeq as defined in the question, DSolve indeed gives an incorrect answer,

sol = DSolveValue[myeq == 0, f0, z];
FullSimplify[(myeq == 0) /. f0 -> %]
(* ((-1 + z) (z^2 - 5 (-2 + z) z Log[1 - z] + (14 + z (-14 + 3 z)) Log[1 - z]^2))/z == 0 *)


which is not, in general, True, as can be seen by evaluating the expression for a few values of z. Note that none of the constants of integration appear in this expression, indicating that the error is not associated with the homogeneous equation. This also can be demonstrated by solving the homogeneous equation directly.

DSolveValue[Hom == 0, f0, z]
FullSimplify[(Hom == 0) /. f0 -> %]
(* Function[{z}, C[1]/z^4 - (C[2] Log[1 - z])/z^4 +
(C[3] (1/(1 - z) - z - 1/4 (-1 - 2 Log[-1 + z])^2))/z^4] *)
(* True *)


Somewhat surprisingly (to me, at least), the correct answer can be obtained as follows. Restructure InHom as

InHomsim = Collect[InHom, _Log, Simplify]
(* -((24 z (24 + z (-48 + z (23 + 2 z))))/(-1 + z)) +
48 (18 + z (-24 + 5 z)) Log[1 - z] + (144 (-2 + z) (-1 + z) Log[1 - z]^2)/z *)


Next, solve

sol1 = DSolveValue[Hom + InHomsim[[1]] == 0, f0, z];
FullSimplify[(Hom + InHomsim[[1]] == 0) /. f0 -> %]
(* True *)
sol2 = DSolveValue[Hom + InHomsim[[2]] == 0, f0, z];
FullSimplify[(Hom + InHomsim[[2]] == 0) /. f0 -> %]
(* True *)
sol3 = DSolveValue[Hom + InHomsim[[3]] == 0, f0, z];
FullSimplify[(Hom + InHomsim[[3]] == 0) /. f0 -> %]
(* True *)


Thus, correct solutions for each of the three additive components of InHomsim can be obtained without difficulty. Simply add them to obtain the correct answer to the original equation.

Collect[(sol1[[2]] + sol2[[2]] + sol3[[2]]) /. C[i_] -> C[i]/3, C[_], FullSimplify];
soltr = Function[{z}, Evaluate@%]
FullSimplify[(myeq == 0) /. f0 -> soltr]
(* Function[{z}, C[1]/z^4 - (C[2] Log[1 - z])/z^4 +
(C[3] (1/(1 - z) - z - 1/4 (1 + 2 Log[-1 + z])^2))/z^4 +
(1/((-1 + z) z^4)) 4 (-6 (7 + 6 \[Pi]^2 (-1 + z) - z (7 + 4 z)) -
2 (-1 + z) Log[1 - z]^3 - 6 Log[1 - z] (-16 + z (15 + 2 z) +
(-1 + z) (-23 + Log[-1 + z]) Log[-1 + z]) + (-1 + z) Log[-1 + z]
(12 + Log[-1 + z] (-69 + 4 Log[-1 + z])) - 6 (-1 + z) Log[1 - z]^2 (5 + Log[z]) +
12 (-1 + z) Log[1 - z] PolyLog[2, z] - 12 (-1 + z) PolyLog[3, 1 - z])] *)
(* True *)


Evidently, DSolve can handle each of the three components of InHomsim but not the whole expression at once.

Why the error?

DSolve also can solve the general equation, Hom + g[z] == 0, correctly.

solgen = DSolveValue[Hom + g[z] == 0, f0, z];
FullSimplify[(Hom + g[z] == 0) /. f0 -> %]
(* True *)


Not surprisingly solgen contains three integrals over g[z] , two of them fairly complicated. Because the g[z] in the question (namely InHom) and (as it happens) the kernels of the integrals both contain expressions with branch cuts, DSolve needs to select the proper contours along which to perform the integrals. Probably, DSolve chose incorrect contours.

• thanks, it does work now. I thought that Assuming would take care of the branch cuts but it's not the case. Only one thing, you have a sol1 written instead of a sol3 in your answer. Thanks again! Commented Nov 27, 2019 at 21:55