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?