We can use the formula
$$F(x)=e^{λx} [ ρ-λμ-\dfrac{1}{2} λ^2 σ^2 ]^{-1}. (1) $$
to derive an expression for F(x) when f(x) is any integer power $x^n$. Begin by observing that for the exponential flow function,
$$F(x)=E[\int_0^\infty\sum_0^\infty\dfrac{1}{n!} (λx_t )^n e^{-ρt} dt|x_0=x]$$ $$F(x)=\sum_0^\infty\dfrac{λ^{n}}{n!} E[\int_0^\infty x_t^n e^{-ρt} dt|x_0=x]$$
Now we can expand the right hand side of (1) in powers of $\lambda$. Since the resulting power series must equal the series just above for all $\lambda$ in a neighborhood that includes $\lambda=0$, we can equate the coefficients of $\lambda^{n}$ on the two sides. This yields an expression for the expected present value of the nth power of x.
To expand the right hand side of (1), note that
$$e^{λx}= \sum_0^\infty\dfrac{λ^{n}}{n!} λ^n x^n,$$
and
$$[ρ-μλ-1/2 σ^2 λ^2 ]^{-1}=\dfrac{1}{\rho}\sum_0^\inftyρ^{-n} (μλ+1/2 σ^2 λ^2 )^n ,$$
For each n, the nth power of the sum of terms in $\lambda$ and $\lambda^2$ in the second of these series must be written out using the binomial theorem. Then the two series can be multiplied together and the coefficients of like powers of $\lambda$ collected to express the product as a single power series. But the first few powers are not too hard, and yield the following results for the discounted present values of x and $x^2$:
$$ E\int_0^\infty x_t e^{-ρt} dt=\dfrac{μ}{ρ^2} +\dfrac{x}{ρ}, (2)$$
$$ E\int_0^\infty x_t^2 e^{-ρt} dt=\dfrac{σ^2}{ρ^2} +\dfrac{2μ^2}{ρ^3}+\dfrac{2μx}{ρ^2} +\dfrac{x^2}{ρ}. (3)$$
Did i not understand how equation (2)and (3), they were obtained?