Which is the covariance matrix between X ~ exponential (lambda) and X squared? I need to find the covariance matrix between $X \sim exp(\lambda)$ and its square $X^2$.
Can I interpret $X^n$ as an Erlang distribution (for which I can have the expectation/variance on Wikipedia :-) ) with parameters $\lambda$ and $n$ and hence
$var(X^2) =  2/\lambda^2$ and $cov(X,X^2) = E[X^3] - E[X]*E[X^2] = \frac{3}{\lambda} - \frac{1}{\lambda} * \frac{2}{\lambda} = \frac{3\lambda -2}{\lambda^2}$ ?

I don't think you need to use the Erlang distribution for that.
As you said: $$cov(X, X^2) = E(X\times X^2) - E(X)E(X^2) .$$
Now, if $X \sim \mathcal{E}xp(\lambda)$, you have that: $$E(X^k) = \frac{k!}{\lambda^{k}}$$ (this is easily shown using the moment generating function of $X$).
So in the end:
$$cov(X, X^2) = \frac{6}{\lambda^3} - \frac{1}{\lambda}\frac{2}{\lambda^2} = \frac{4}{\lambda ^3}$$
