Ok here's a possible linear transformation way, I think it'll work:
So I'll just write the basis in two different ways:
B={\(e^x,xe^x,x^2e^x,x^3e^x,\dots \)} = {\(y_0,y_1,y_2,y_3, \dots\)}
\[Dy_n = y_n+ny_{n-1}\]
So you can try that out, and then see we have two matrices sorta working here, the first is an identity matrix and the second is sorta like a shift matrix and has n in it. This is actually just the regular derivative matrix for if our basis didn't have \(e^x\) in there (you can kinda believe this from the product rule just multiplying \(e^x\) and then differentiating the polynomial):
\[Dy_n = I y_n+Ay_n\]
So to be clear about what I'm writing:
\[Iy_n = y_n\]\[Ay_n = n y_{n-1}\]
\[D=I+A\]
Ok so the problem we have is that \(A\) is an infinite matrix. Luckily though, we don't need to know all the entries, just enough up to what we're interested in, since we can see that the derivative matrix will only map from our basis at the place \(n\) to \(n-1\) that means the inverse matrix which exactly reverses this can at most map backwards to \(n\) and \(n+1\). So that means we just need to find the derivative matrix up to \(n+1\) more than what we want to integrate, simple as that.
I'll show an example.
\[\int x^2e^x dx\]
So we have to calculate the derivative matrix for up to 3. So that means we have:
\[D = I+ A = \begin{bmatrix}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{bmatrix} + \begin{bmatrix}
0& 1 & 0 & 0\\
0 & 0& 2& 0\\
0 & 0 & 0 & 3\\
0 & 0 & 0 & 0
\end{bmatrix} =\begin{bmatrix}
1& 1 & 0 & 0\\
0 & 1& 2& 0\\
0 & 0 & 1 & 3\\
0 & 0 & 0 & 1
\end{bmatrix} \]
Now we invert this matrix... with matlab...
\[ D^{-1} = \begin{bmatrix}
1& -1 &2 & -6 \\
0 & 1& -2& 6\\
0 & 0 & 1 & -3\\
0 & 0 & 0 & 1
\end{bmatrix} \]
Multiply by the vector representing our integral:
\[\begin{bmatrix}
1& -1 &2 & -6 \\
0 & 1& -2& 6\\
0 & 0 & 1 & -3\\
0 & 0 & 0 & 1
\end{bmatrix} \begin{bmatrix}
0\\
0\\
1\\
0
\end{bmatrix} = \begin{bmatrix}
2\\
-2\\
1\\
0
\end{bmatrix}\]
Ok moment of truth, hopefully I didn't mess this up...
\[\int x^2e^xdx = e^x(2-2x+x^2) + C\]
Hey it checks out :P
Also, maybe we can possibly find a general form for the inverse of this matrix?