I think that the substitution $x=\tan(a)$ is from far away the simplest to use for the calculation of $$I_m=\int \frac{\mathrm{d}x}{(x^2+1)^{m}}=\int \cos^{2(m-1)}(a)\,\mathrm{d}a$$ At this point, the only thing to be done is a power reduction of the cosine which will transform the cosine raised to power into a linear combination of cosines of multiple angles; this is what you did.
The general formulas are quite simple. If $m$ is odd,
$$\cos^m(a)=\frac 1{2^{m-1}}\sum_{k=0}^{\frac{m-1}2}\binom{m}{k} \cos \Big( (m-2 k)\,a\Big)$$ and, if $m$ is even,
$$\cos^m(a)=\frac 1{2^{m}}\binom{m}{\frac m2} +\frac 1{2^{m-1}}\sum_{k=0}^{\frac{m}2-1}\binom{m}{k} \cos \Big( (m-2 k)\,a\Big)$$ and the integration does not make any problem.
You could fine similar reduction formulas for $\sin^m(a)$.
If you do not use this substitution, you will face hypergeometric functions (as Dmoreno already mentioned) and $$I_m=x \, _2F_1\left(\frac{1}{2},m;\frac{3}{2};-x^2\right)$$ which can be represented as $$I_m=\frac{P_{2m-3}(x)}{(1+x^2)^{m-1}}-k_m \tan^{-1}(x)$$ where $P_n(x)$ represents a polynomial of degree $n$.