How to compute the indefinite integrate of n-th tetration of x? How can the following indefinite integral be computed ?
$ \int{x↑↑n} dx $
where $n$ = {$x$ $\in$ $N^+$ : ${x > 2}$}
Here ${x↑↑n}$ refers to $n$th tetration of $x$. I tried searching over the internet to find a general formulae for calculating it when the value of  $n$ varies and came across a single paper which had an application of W-Lambart function. 
Is there any other method to compute it ?
 A: The solution expressed on the form of series isn't generally accepted as a closed form.  However, the choice of what to call closed-form and what not is rather arbitrary. In common sens, a closed form is a combination of a $\underline{\text{finite}}$ number of $\underline{\text{standard}}$ functions. 
Then the question is what is a standard function ? When a function (defined by an infinite series or by an integral or by other means) becomes a standard function ? 
For an example related to $ \int{x↑↑n} dx $ in the paper : https://fr.scribd.com/doc/34977341/Sophomore-s-Dream-Function the function named Sphd covers a lot of cases.
$$\text{Sphd}(\alpha,x)=\int_0^x t^{\alpha t}dt$$
generalized with :

$$ \int_0^x{x↑↑n} \: dx = \text{Sphd}_{n}(1,...;x)\qquad \text{Eq.}(12:5)$$
And from Eq.$(12:6)$
$$ \int_a^b{x↑↑\infty} \: dx =\text{Sphd}_{\infty}(1,...;b)-\text{Sphd}_{\infty}(1,...;a)=-\int_a^b\frac{\text{W}(-\ln(t))}{\ln(t)}dt$$
W is the Lambert W function.
Of course, all this is purely formal since Sphd isn't a standardized special function. 
Creating symbols for new functions may appear as vicious circle. A discussion about the definition of new special functions, standardization of such functions and further use can be found in :  https://fr.scribd.com/doc/14623310/Safari-on-the-country-of-the-Special-Functions-Safari-au-pays-des-fonctions-speciales . 
A: Your can try to compute it approximately.
As an example. If you have the integral.
$$
I_2(a)=\int_1^a x^{x^x} dx\,,
$$ with $a\gg1$. In fact this faction grows very rapidly, it will be works for a=3, with enough good accuracy.
We can represent this integral in the following form
$
I_2(a)=\int_1^a e^{S(x)} dx
$
and do the integration by parts
$$
I(a)=\int_1^a e^{S(x)} dx=\int_1^a \frac{de^{S(x)}}{S'(x)} dx=\frac{S(a)}{S'(a)}-\frac{S(1)}{S'(1)}-\int_1^a dxe^{S(x)} \frac{d}{dx} \frac{1}{S'(x)}. 
$$
The main contribution is given by the term $\frac{S(a)}{S'(a)}$, the integral term is a the next correction. You can continue to do the integration by part to improve your accuracy.
In such case 
$I_2(a)=\dfrac{a^{a^a}}{a^a(1/a+\ln a(\ln a+1 ))}$
The I_2(3)=$1,07 *10^{11}$, the result by numerical integration is $1,11*10^{11}$.
In fact this asymptotics depends only on upper limit of integration.
In the general case, using the notation presented above we have
$$
I_n(a)=\int_0^a{x↑↑n} dx=\frac{x↑↑n}{d/dx({x↑↑(n-1)} \ln(x))}|_{x=a}
$$
