# DSolve apparently returns wrong symbolic result

Bug introduced in 9.0 and persisting through 11.0.1 or later

I am trying to get Mathematica to calculate the correct solution $x(t)$ of the following ODE:

$\dot x = a x(t)+ct^n$

for $a>0$, $n\ge 1$, $c>0$ and $x(0) = x_0$.

My final goal is to calculate the asymptotic behavior for $t\rightarrow \infty$ (i.e. the fastest growing term). I can do so on a piece of paper and I get

$\left(x_0+\frac{c n!}{a^{n+1}}\right)e^{at}$.

I am trying to do the same thing using Mathematica with this code:

eqn = D[x[t], t] == a*x[t] + c*t^n;
sol = DSolve[{eqn, x[0] == x0}, x, t]
Simplify[eqn /. sol] (*1*)
solX[t_] = (x[t] /. sol[[1]]) (*2*)
asymptX[t_] =
Limit[solX[t]/(Exp[a t]), t -> \[Infinity],
Assumptions -> {a > 0, n >= 1, c > 0}]*Exp[a t] (*3*)


The line marked "(*1*)" returns {True}, so Mathematica thinks its solution is correct. The full solution (*2*) is given as

$\frac{e^{a t} (a t)^{-n} \left(-c t^n \Gamma (n+1,a t)+a c \Gamma (n+1,0) (a t)^n+a x_0 (a t)^n\right)}{a}$

and the asymptotic solution (*3*) that is returned is

$e^{a t} (c \Gamma (n+1)+x_0) = e^{a t}(cn! + x_0)$.

As you can see, the problem is that the factor $1/a^{n+1}$ is missing from the asymptotic solution compared to my by-hand calculation. I also calculated the numerical solution with Mathematica and then checked and compared all the solutions:

asymptXReal[t_]=(x0+c n!/a^(n+1))*Exp[a t]
Block[{a = 0.8,n=5,c=0.123,x0=0.123},
nsol=NDSolve[{eqn,x[0]==x0},x,{t,0,15}];
Plot[{Evaluate[x[t]/.nsol],solX[t],asymptXReal[t],asymptX[t]},
{t,0,15},PlotLegends->Placed["Expressions",Below]]]


which yields the following graph:

You see that my manual derivation agrees asymptotically with the numerical solution, while Mathematica's full and asymptotic symbolic solutions are both off (but consistent with each other). Apparently, the asymptotic solution is wrong because the full solution is wrong to begin with. But why is that? I guess I am just making a stupid mistake here. Can anyboday help me find it?

• @xzczd, thanks for the diligence in version tracking! :) – J. M. will be back soon Jun 28 '16 at 10:23

This turns out to be a bug since v9. v8.0.4 gives the correct result:

\$Version
eqn = D[x[t], t] == a*x[t] + c*t^n;
bc = x[0] == x0;
sol = Simplify[DSolve[{eqn, bc}, x[t], t], {a > 0, n >= 1, c > 0}]
(* "8.0 for Microsoft Windows (64-bit) (October 24, 2011)" *)
(* {{x[t] -> a^(-1 - n) E^(a t) (a^(1 + n) x0 + c Gamma[1 + n] - c Gamma[1 + n, a t])}} *)


The analytic solution given in v9 and v10 isn't in consistent with the b.c.:

Simplify[x[0] == x0 /. sol, {a > 0, n >= 1, c > 0}]
(* {a == 1} *)


If you're in v10, a possible workaround is to add Assumptions to DSolve:

sol2 = DSolve[{eqn, x[0] == x0}, x, t, Assumptions-> {a > 0, n >= 1, c > 0}]
(* {{x -> Function[{t},
a^(-1 - n) E^(a t) (a^(1 + n) x0 + c Gamma[1 + n] - c Gamma[1 + n, a t])]}} *)


Another workaround that also works in v9 is to find the general solution first and then solve for the constant:

sol3 = # /. Solve[#[[1, 1, -1]][0] == x0, C[1]][[1]] &@DSolve[eqn, x, t];

Simplify[sol3 /. Function -> func, {a > 0, n >= 1, c > 0}] /. func -> Function
(*{{x -> Function[{t}, E^(a t) (x0 + c Gamma[1 + n] - a^(-1 - n) c Gamma[1 + n, a t])]}}*)

• Great idea to check the b.c.! It was kind of mysterious to me how the bug could affect DSolve and Simplify at the same time to make Mathematica think its solution was correct. Thanks for trying the different versions and the workaround. I have sent a bug report to Wolfram support. – Oguk Jun 28 '16 at 7:06
• @Oguk What was the outcome of your bug report? – bbgodfrey Aug 13 '16 at 3:53