Evaluation of $\displaystyle \int x^8\sqrt{x^2+1}\mathrm dx$
$\bf{My\: Try:}$ Put $x=\tan \theta\;,$ Then $\mathrm dx = \sec^2 \theta\mathrm d \theta$
So $$I = \int \sec ^{3}\theta \cdot \tan^{8}\theta\mathrm d \theta=\int \frac{\sin^{8}\theta}{\cos^{5}\theta }\mathrm d\theta$$
Now, how can I evaluate this integral? Please help.
Apply the reduction formula $$\int x^n\sqrt{1+x^2}dx= I_n=\frac{x^{n-1}}{n+2}(1+x^2)^{3/2}-\frac{n-1}{n+2}I_{n-2} $$ 4 times to reduce the integral $I_n$ to the manageable $I_0$
\begin{align} \int x^8\sqrt{1+x^2}\ dx =7\left(\frac{x^7}{70}-\frac{x^5 }{80}+\frac{x^3 }{96}-\frac{x}{128}\right)(1+x^2)^{3/2}+\frac7{128}I_0 \end{align} where $ I_0=\int \sqrt{1+x^2}\ dx =\frac{x}{2}(1+x^2)^{1/2}+\frac{1}{2}\sinh^{-1}x $
By substituting $x=\sinh\theta$
$$ I = \int \cosh(\theta)^2 \sinh(\theta)^8\,d\theta = \int\sinh(\theta)^8\,d\theta+\int\sinh(\theta)^{10}\,d\theta$$ and the last integrals are straightforward to compute through De Moivre's formula, since
$$ J(n) = \int \sinh^{2n}(\theta)\,d\theta = \frac{1}{4^n}\sum_{k=0}^{2n}\binom{2n}{k}\int (-1)^k e^{-k\theta}e^{(2n-k)\theta}\,d\theta $$ gives: $$ J(n) = \frac{(-1)^n \theta}{4^n}\binom{2n}{n}+\frac{2}{4^n}\sum_{k=0}^{n-1}\binom{2n}{k}\frac{(-1)^k \cosh((2n-2k)\theta)}{2n-2k}.$$
Hint: Using the identity $\tan^2\theta=\sec^2\theta-1$ you can rewrite your integral as $$\int{\sec^3(\theta)(\sec^2{(\theta)}-1)^4}\ d\theta$$
HINT:
Integrating by parts,
for $m+1\ne0,$
$$I_m=\int x^m\sqrt{x^2+1}\ dx=\sqrt{x^2+1}\int x^m\ dx-\int\left(\dfrac{d\sqrt{x^2+1}}{dx}\int x^m\ dx\right)dx$$
$$(m+1)I_m=x^{m+1}\sqrt{x^2+1}-\int\dfrac{x^{m+2}}{\sqrt{x^2+1}}dx$$
$$=x^{m+1}\sqrt{x^2+1}-\int\dfrac{x^m(1+x^2)-x^m}{\sqrt{x^2+1}}dx$$
$$(m+2)I_m=x^{m+1}\sqrt{x^2+1}+I_{m-2}-\int\dfrac{x^{m-2}}{\sqrt{x^2+1}}dx$$
The last part can be reduced further as $x^{m-2}=x^{m-4}(1+x^2)-x^{m-4}$
You can use Euler substitution: $\sqrt{x^2+1}=-x+t$. Then you get the integral:
$\int \sqrt{x^2+1}x^8dx = \int t*(\dfrac{t^2-1}{2t})^8*\dfrac{t^2+1}{2t^2}dt$
$\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{\mathrm{i}} \newcommand{\iff}{\Leftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\, #2 \,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$
\begin{align} \color{#f00}{\int\root{x^{2} + 1}x^{8}\,\dd x} & = {1 \over 9}\int\root{x^{2} + 1}\,\dd\pars{x^{9}} \\ & = {1 \over 9}\,x^{9}\root{x^{2} + 1} - {1 \over 9}\int{x^{10} \over \root{x^{2} + 1}}\,\dd x \end{align}
With the sub$\ldots\quad$ $\ds{t \equiv \root{x^{2} + 1} - x\quad\imp\quad x = {1 - t^{2} \over 2t}}$: \begin{align} & \color{#f00}{\int\root{x^{2} + 1}x^{8}\,\dd x} = {1 \over 9}\,x^{9}\root{x^{2} + 1} - {1 \over 9216}\int{\pars{1 - t^{2}}^{10} \over t^{11}}\,\dd t \\[3mm] & = {1 \over 9}\,x^{9}\root{x^{2} + 1} - {1 \over 9216}\sum_{n = 0 \atop {\vphantom{\large A}n \not= 5}}^{10}{10 \choose n}\pars{-1}^{n} \int t^{2n - 11}\,\dd t + {7 \over 256}\ln\pars{t} \\[3mm] & = {1 \over 9}\,x^{9}\root{x^{2} + 1} - {1 \over18432}\sum_{n = 0 \atop {\vphantom{\large A}n \not= 5}}^{10} {10 \choose n}{\pars{-1}^{n} \over n - 5}\,t^{2n - 10} + {7 \over 256}\ln\pars{t} \\[3mm] & = \color{#f00}{{1 \over 9}\,x^{9}\root{x^{2} + 1} - {1 \over18432}\sum_{n = 0 \atop {\vphantom{\large A}n \not= 5}}^{10} {10 \choose n}{\pars{-1}^{n} \over n - 5}\,\pars{\root{x^{2} + 1} - x}^{2n - 10}} \\[3mm] & \color{#f00}{\phantom{==}+ {7 \over 256}\ln\pars{\root{x^{2} + 1} - x}} \end{align}
Expanding my comment:
setting $x=\sinh t$, $dx=\cosh t\ dt$, the integral becomes $$\int\sqrt{1+\sinh^2 t}\sinh^8 t\cosh t\ dt=\int\sinh^8t\cosh^2t\overset{P.I.}{=}\frac19\sinh^9t\cosh t-\int \sinh^{10}t\ dt$$ To evaluate this last integral you can use the definition $$\sinh t=\frac{e^t-e^{-t}}2$$ so that $$\sinh^{10} t=\frac{1}{2^{10}}\sum_{k=1}^{10}{10 \choose k}(-1)^ke^{(10-k)t}e^{-kt}=\frac{1}{2^{10}}\sum_{k=1}^{10}{10 \choose k}(-1)^ke^{(10-2k)t}$$
but I leave this up to you
Remark: of course you can substitute $\cosh^2 t$ in the second integral too, so that you then evaluate two $\sin^m t$ integrals. As done in the answer below
Let $x=\sinh t.$ Then $\sqrt {1+x^2}=\cosh t.$ And the integral is $$\int \cosh^2 t \sinh^8 t \;dt=\int (1+\sinh^2 t)\sinh^8 t dt=\int (\sinh^8 t+\sinh^{10}t) \;dt.$$ To evaluate $\int \sinh^{2 n}t \;dt $: For brevity let $S=\sinh t$ and $C=\cosh t,$ and let $\int \sinh^{2 n}t \;dt =V(2 n).$
Then $$V(2 n)=\int S^{2 n-1}\;dC=$$ $$=S^{2 n-1}C-\int C\;d(S^{2 n-1})=$$ $$=S^{2 n-1}C-\int (2 n-1)C^2 S^{2 n-2}\;dt=$$ $$=S^{2 n-1}C-\int (2 n-1)(1+S^2)S^{2 n-2}\;dt=$$ $$=S^{2 n-1}C-(2 n-1)(V(2 n-2)+V(2 n)).$$ This implies $$2 n V(2 n)=S^{2 n-1}C-(2 n-1)V(2 n-2).$$
Alternative Euler approach:
$$\begin{align*} I &= \int x^8 \sqrt{x^2+1} \, dx \\ &= \frac19 x^9 \sqrt{x^2+1} - \frac19 \int \frac{x^{10}}{\sqrt{x^2+1}} \, dx \tag1\\ &=\frac19 x^9 \sqrt{x^2+1} - \frac1{18} \int \frac{y^5}{\sqrt{y^2+y}} \, dy \tag2\\ &=\frac19 x^9 \sqrt{x^2+1} - \frac19 \int \frac{z^9}{(1-2z)^6} \,dz \tag3\\ &=\frac19 x^9 \sqrt{x^2+1} + \frac1{9216} \int \frac{(1-w)^9}{w^6} \,dw \tag4\\ \end{align*}$$
and proceed with partial fraction expansion.
The given integral falls into case $4$ of Chebyshev's differential binomial, hence it can be evaluated by substituting $t=\sqrt{1+\frac1{x^2}}$. Then, the integral becomes
$$\int\frac{\mathrm dt}{(t^2-1)^6}-\int\frac{\mathrm dt}{(t^2-1)^5}$$
which can be finished off by PFD.