First, I'll generalize the problem a bit by taking $A_{jk}=\int_{\mathbb{R}} x^{j+k}\,d\mu$ where $\int_{\mathbb{R}}d\mu=1$
(formally, $\mu$ is a probability measure and $d\mu/dx$ is the probability density function. but this is basically just notation so that I don't have to worry about which particular inner product I'm doing.)
Since I can consider this as $A_{jk}=\langle x^j,x^k\rangle$ where $\langle f(x),g(x)\rangle =\int_\mathbb{R} f(x) g(x)\,d\mu$, I can do Gram-Schmidt on this inner product to get a family of orthogonal polynomials $\{P_k(x)\}_{k=0,1,2,\cdots}$ where $\text{deg }P_k=k$ and $\langle P_j(x),P_k(x)\rangle =\delta_{jk}$ i.e. orthonormal
(actually, I didn't need to require $\int_{\mathbb{R}}d\mu=1$; all I should have said was that $\mu(x)$ be non-decreasing on its support and that $\int_\mathbb{R}f(x)d\mu$ be finite for any polynomial $f(x)$ )
As an aside, the case at hand corresponds to $d\mu = \chi_{(0,1)}\,dx$ i.e. uniform support on $[0,1]$. This corresponds to the shifted Legendre polynomials; by contrast, the case $d\mu = \chi_{(-1,1)}\,dx$ corresponds to regular Legendre polynomials.
Anyways. Once we have the family of OP's, we can formally expand any monomial $x^j$ in this family as $x^j = \sum_n C_{jn} P_n(x)$. Using the orthonormality of the inner product, we further have $\langle x^j,P_k\rangle = \sum_n C_{jn} \langle P_n , P_k\rangle = C_{jk}$
Hence computing the matrix $C_{jk}$ amounts to computing inner products of the monomials with the OPs.
The reason this is a smart idea is because now one can use this to express the original matrix $A$:
$$A_{jk}=\langle x^j,x^k \rangle = \sum_{mn}C_{jm}C_{kn}\langle P_m,P_n\rangle = \sum_{mn}C_{jm}C_{kn}\delta_{mn} = \sum_{n} C_{jn}C_{kn}.$$
If I regard $C$ as a matrix with elements $C_{jk}$, though, I can interpret this as matrix multiplication. Hence $A=CC^T$.
Additionally, one should be able to write any monomial $x^j$ using only polynomials of degree $k\leq j$. Hence $C_{jn}=\langle x^j,P_n = 0$ if $n>j$, so that $C$ is a lower-triangular matrix.
Hence $A$ can be explicitly written as an $LU$-factorization with $U=L^T$.
If we further require that the orthogonal polynomials be monic i.e. $P_k(x)=x^k+\cdots$ then we have $C_{kk}=\langle P_k,x^k \rangle =1$
Hence $C$ can be written as $C=I+M$ where $M$ is strictly lower-triangular
The only task now is to find $C^{-1}=(I+M)^{-1}$; if we can do this, then it's easy to see that $A^{-1}=(C^{-1})^T C^{-1}$.
The main computational task is to find $C_{jk}$ in a nice way.
Actually, it looks like I can't do all these assumptions at once. If $P_n(x)$ is an orthonormal sequence, then I can't require it be monic as well.
So probably best to ignore everything past the LU-factorization statement.
