I am not entirely sure where I went wrong with my approach so I wanted to share it with you guys, thanks for any insights or for spotting some silly mistake.
$$\int \frac{dx}{\tan(ax)^n}$$
First, I made a substitution as such: $u = \tan(ax)$ thus obtaining $$a\sec(ax)^2dx = du$$ rewriting it gives us $$a(\tan(ax)^2+1)dx = du$$
Then I rewrote the integral with the substitution as follows:
$$\frac{1}{a}\int \frac{du}{u^n(u^2 + 1)}$$
Then I integrate as follows (following step is wrong, this should be done with partial decomposition but is it possible?):
$$\frac{1}{a}\int \frac{du}{u^n(u^2 + 1)} = \frac{1}{a}\ln \left| u^{n+2} + u^n \right| \cdot \frac{1}{(n+1)u^{n+1} + (n-1)u^{n-1}} + C$$
Then substituting back results in:
$$\frac{1}{a}\ln \left| u^{n+2} + u^n \right| \cdot \frac{1}{(n+1)u^{n+1}} + (n-1)u^{n-1} + C = \frac{1}{a}\ln \left| \tan(ax)^{n+2} + \tan(ax)^n \right| \cdot \frac{1}{(n+1)\tan(ax)^{n+1} + (n-1)\tan(ax)^{n-1}} + C$$
Testing my answer with $\cot(3x)^4$ I get a very similar graph to the answer in the book but it is slightly off I can't seem to find the mistake. Thanks for any insights.
It seems that my mistake is the integration step as I do not get a log integral, it is wrong. I should consider partial decomposition, but then is it possible to generalise it?
$\textbf{Update:}$ so far I have only realised that if we play around with different powers $n$ we notice that if there is a general expression, it is different for odd and even powers; I conjecture that we indeed can find it, but two separate expressions.
$$I_n=\int \frac{du}{u^n(u^2 + 1)}$$
If you use Mathematica, it returns, as @138Aspen wrote in comments $$I_n=\frac{u^{1-n} }{1-n}\, _2F_1\left(1,\frac{1-n}{2};\frac{3-n}{2};-u^2\right)$$ or $$I_n=-\frac{1}{2} i^{n+1}\, B_{-u^2}\left(\frac{1}{2}-\frac{n}{2},0\right)$$
The problem is that, if $n$ is odd, Mathematica returns the unpleasant
ComplexInfinity.I think that it would be better to analyze the structure of the individual integrals $$\large\color{blue}{I_{2n}=(-1)^n \tan^{-1}(u)+P_n(u)}$$ where $$\left( \begin{array}{cc} n & P_n \\ 1 & -\frac{1}{u} \\ 2 & -\frac{1}{3 u^3}+\frac{1}{u} \\ 3 & -\frac{1}{5 u^5}+\frac{1}{3 u^3}-\frac{1}{u} \\ 4 & -\frac{1}{7 u^7}+\frac{1}{5 u^5}-\frac{1}{3 u^3}+\frac{1}{u} \\ 5 & -\frac{1}{9 u^9}+\frac{1}{7 u^7}-\frac{1}{5 u^5}+\frac{1}{3 u^3}-\frac{1}{u} \\ 6 & -\frac{1}{11 u^{11}}+\frac{1}{9 u^9}-\frac{1}{7 u^7}+\frac{1}{5 u^5}-\frac{1}{3 u^3}+\frac{1}{u} \\ \end{array} \right)$$
$$\large\color{blue}{I_{2n+1}=(-1)^n \frac{1}{2} \log \left(\frac{u^2+1}{u^2}\right)+Q_n(u)}$$ where $$\left( \begin{array}{cc} n & Q_n \\ 1 & -\frac{1}{2 u^2} \\ 2 & -\frac{1}{4 u^4}+\frac{1}{2 u^2} \\ 3 & -\frac{1}{6 u^6}+\frac{1}{4 u^4}-\frac{1}{2 u^2} \\ 4 & -\frac{1}{8 u^8}+\frac{1}{6 u^6}-\frac{1}{4 u^4}+\frac{1}{2 u^2} \\ 5 & -\frac{1}{10 u^{10}}+\frac{1}{8 u^8}-\frac{1}{6 u^6}+\frac{1}{4 u^4}-\frac{1}{2 u^2} \\ 6 & -\frac{1}{12 u^{12}}+\frac{1}{10 u^{10}}-\frac{1}{8 u^8}+\frac{1}{6 u^6}-\frac{1}{4 u^4}+\frac{1}{2 u^2} \\ \end{array} \right)$$
The patterns are quite obvious.