I did this by using Feynman trick and differential equation. I guess there are a few functions presented as an integral, for which this differentiation trick does not hold true owing to the absence of derivative on some set; or owing to humans error – we all learn by making mistakes. Anyway, take a look at my solution (my attempt at solving the problem I post here!). We define
$$I(n) = \int_0^1 x^{nx}dx$$
Step 1. I take the derivative of the integral with respect to $n$, getting
$$\frac{dI}{dn} = \int_0^1 x^{nx} \cdot \ln x \cdot x \cdot dx$$
Step 2: Next, I notice that $(x^{nx})^\prime = n x^{nx} \left(\ln \left(x\right) +1\right)$ and using this I try to maintain the integral through the differential method:
$$\int_0^1 x^{nx} \ln \left(x\right) \cdot xdx = \frac{1}{n}\int_0^1 n \cdot x^{nx} \left(\ln \left(x\right) + 1 – 1 \right)xdx $$
So, I get
$$ \frac{dI}{dn} = \frac{1}{n} \int_0^1 n \cdot x^{nx} \left(\ln \left(x\right) + 1\right) xdx – \int_0^1 x^{nx} xdx$$
Note that I reduce $1/n \cdot n = 1$ in the second integral while in the other integral I keep it.
Step 3: I integrate the first integral by parts (you discerned my preparation for this during the previous step!); for that, I put $u = x$, then $du = dx$; and $dv = n\cdot x^{nx}\left(\ln \left(x\right) + 1\right)dx$, then $v = x^{nx}$. So, I have
$$\int_0^1 n \cdot x^{nx} \left(\ln \left(x\right) + 1\right) xdx = x \cdot x^{nx} |_0^1 – \int_0^1 x^{nx} dx = 1 – I(n)$$
All steps before the following ones are considered (by me) to be correct and have thoroughly been verified; if my solution is incorrect, then I failed to solve my problem during the next problems!
Step 4: The second integral should also be integrated by parts; for that, I put $u = x$, then $du = dx$; and $dv = x^{nx}dx$. As an antiderivative of $x^{nx}$, I take $v = \int_0^x t^{nt}dt$. Thus,
$$\int_0^1 x^{nx} xdx = \left(x \int_0^x t^{nt} dt\right)|_0^1 – \int_0^1 \left(\int_0^x t^{nt}dt\right)dx$$
Note that
$$\left(x \int_0^x t^{nt} dt\right)|_0^1 = 1 \cdot \int_0^1 t^{nt}dt – 0 \cdot \int_0^0 t^{nt}dt = I(n)$$
Next, we try to switch the order of integration:
$$ \int_0^1 \int_0^x t^{nt} dt dx = \int_0^1 t^{nt}dt \int_0^t dx = \int_0^1 t^{nt} tdt$$
Then
$$ \int_0^1 x^{nx} x dx = \frac{1}{2} I(n)$$
Subsequently, the initial integral should satisfy the equation
$$ \frac{dI}{dn} = \frac{1}{n} – \frac{I}{n} – \frac{I}{2} = \frac{1}{n} – \frac{n+2}{2n}I(n)$$
Or, equivalently,
$$\frac{dI}{dn} + \frac{n + 2}{2n}I = \frac{1}{n}$$
which is the linear homogenous differential equation of the first order; in order to solve this, I used the method called Lagrange: in accordance with it, I had to solve the corresponding inhomogeneous equation $I^\prime + \varphi(n)I = 0$ and then vary the constant $C = C(n)$; let us dive in details: I had to solve the equation
$$ \frac{dI}{dn} + \frac{n+2}{2n}I = 0$$
or
$$\frac{dI}{I} + \frac{n+2}{2n}dn = 0$$
supposing that neither $I \ne 0$ and $n \ne 0$. Integrating yields:
$$ \ln I + \frac{n}{2} + \ln n = \ln C \Rightarrow I(n) = \frac{C}{n}e^{-\frac{n}{2}}$$
where
$$\int \frac{n + 2}{2n}dn = \int \left(\frac{1}{2}+\frac{1}{n}\right)dn =\frac{1}{n} + \ln n$$
See that $I = 0$ is a solution that might be produced from this one, so it is not an isolated solution and we did not lose it. Next, I vary $C = C(n)$. Then I find the derivative to be
$$\frac{dI}{dn} = \frac{d}{dn} \left(\frac{C(n)}{n}e^{-\frac{n}{2}}\right) = \frac{dC}{dn} \cdot \frac{1}{n} \cdot e^{-\frac{n}{2}} – C(n) \cdot \frac{n+2}{2n^2}e^{-\frac{n}{2}}$$
Putting $I$ and $dI/dn$ into the initial equation, I arrive at another
$$ \frac{dC}{dn} \cdot \frac{1}{n} \cdot e^{-\frac{n}{2}} – C(n) \cdot \frac{n+2}{2n}e^{-\frac{n}{2}} + \frac{n+2}{2n} \cdot \frac{C}{n}e^{-\frac{n}{2}} = \frac{1}{n}$$
or, equivalently,
$$ \frac{dC}{dn} \cdot \frac{1}{n} \cdot e^{-\frac{n}{2}} = \frac{1}{n} $$
From this,
$$C(n) =\int e^{\frac{n}{2}}dn = 2e^{\frac{n}{2}}+A$$
where A is another constant, this time, real constant. In our case, this constant is the specific value for which to be defined, I see that $I(1) = \int_0^1 x^x dx = \alpha$, then from
$$ I(n) = \frac{2}{n} + \frac{A}{n}e^{-\frac{n}{2}}$$
setting $n = 1$, I compute $A = \left(\alpha – 2\right)\sqrt{e}$
The terminal result is
$$ I(n) = \frac{2}{n} + \frac{\left(\alpha – 2\right)\sqrt{e}}{n}e^{-\frac{n}{2}} = \frac{2\sqrt{e^{n-1}}+\alpha -2}{n\sqrt{e^{n-1}}}$$
Besides which, $I(+\infty) = 0$. When $n = 2$ my formula gives the result which is approximately equal to that given by WolframAlpha programme; the difference is small but the error increases when $n$ becomes large, though the integral is still approaching zero at infinity and both formulae coincide.
Links:
The true calculation by Wolfram: https://www.wolframalpha.com/input/?i=%5Cint_0%5E1+x%5E%282x%29dx
Best Answer
You messed up at this step: $$\int_0^1 \int_0^x t^{nt} dt dx = \int_0^1 t^{nt}dt \int_0^t dx = \int_0^1 t^{nt} tdt$$
Iterated integrals where the inner bounds depend on the outer ones can't be decoupled like that. Instead, it should be $$\int_0^1 \int_0^x t^{nt} dt dx = \int_{0}^{1}\int_{t}^{1}t^{nt}dxdt = \int_{0}^{1}t^{nt}\left(1-t\right)dt$$
which doesn't provide any information, as it just makes the integral "loop back" to $$\int_0^1t^{nt}tdt = \int_0^1t^{nt}tdt$$
The statement $$\int_0^1t^{nt}tdt = \frac{1}{2}I(n)$$ and everything after that is not correct then.