Express $I_{n}(x)=\int (x^{2}+1)^{(n-0.5)} dx$ with recursion Let $I_{n}(x)=\int (x^{2}+1)^{(n-0.5)} dx$. I want to express $I_{n}(x)$ with earlier terms $\sum_{k\leq n} h_{k} I_{k}(x)$ where $h_{k}$ is a constant. I have no idea by which terms. 
If I am going to do this manually, I am suspecting that I need to do here some trigonometric tricks to get this formula into some nice and easy form and then just simple integration by parts (here). 
Trial 1. Dead end.
Perhaps, useful $$\int \frac{1}{1+x^{2}} dx = Artcan(x) +C$$. 
Trial 2. Dead end.
Let $II_{n}(x)=\int cos^{n}(x) dx$ so by a few times integration by parts: 
$$n II_{n}(x)-(n-1)II_{n-2}(x)=sin(x)cos^{n-1}(x)$$,
this may become useful -- no idea, look I find dead end here. 
Trial 3 and question. Some computational method?
I would like to know whether some program can help with this kind of problems. I find it very erroneous to do many integration-by-parts again-and-again, some computer software to check calculations or help with this kind of repetitive problems that can be very time-consuming to realize?
 A: Substituting
$$x:=\sinh y\ ,\quad dx=\cosh y\ dy$$
you get
$$I_n(x)=\int\cosh^{2n}y\ dy\Bigr|_{y:={\rm arsinh} x}\ .$$
Now you have to set up a recursion for the integrals
$$J_n(y):=\int\cosh^{2n}y\ dy\ .$$
This can be done using partial integration: 
$$\eqalign{J_n(y)&=\int\cosh^{2n-1}y\ \cosh y\ dy\cr &= \cosh^{2n-1}y\ \sinh y -
(2n-1)\int \cosh^{2n-2}y\ \sinh^2 y\ dy \cr &= \cosh^{2n-1}y\ \sinh y -
(2n-1)\int \cosh^{2n-2}y\ (\cosh^2 y-1)\ dy\cr
&=\cosh^{2n-1}y\ \sinh y -(2n-1)\bigl(J_n(y)-J_{n-1}(y)\bigr)\ .\cr }$$ 
This gives
$$J_n(y)={1\over 2n}\bigl(\cosh^{2n-1}y\ \sinh y +(2n-1) J_{n-1}(y)\bigr)\ .$$
Finally the back-substitution: $\sinh y:=x$, $\ \cosh y:=\sqrt{1+x^2}$.
A: Plugging the integral into Mathematica gives
$$I_n(x) = x\ {}_2F_1\bigl(0.5, 0.5 - n; 1.5; -x^2\bigr)$$
The definition of the hypergeometric function ${}_2 F_1$ is well explained on Wikipedia, including a set of recurrence relations. The one of interest to you, in Wikipedia's notation, is
$$F(b-) = \frac{2b - c + (a-b)z}{b - c} F + \frac{1 - z}{c - b}F(b+)$$
You get this by equating the second and fifth lines of Gauss's contiguous relations in the Wikipedia article. Plugging in your values, if I've done the algebra right, you get
$$I(n) = \frac{0.5 + 2n + nx^2}{1+n}I(n-1) + \frac{1 + x^2}{1 + n}I(n-2)$$
Hypergeometric functions are well studied and many numeric libraries include methods to calculate them reasonably efficiently.
