Just follow the Gram-Schmidt algorithm outlined there, where your inner product is the weighted inner product $$\langle u,v\rangle_w = \int_a^b u(x)v(x)w(x)\,dx.$$ Here, $w(x)$ is called the weight function.
In your context, it sounds like you want to generate the first four Laguerre polynomials, $\{L_0(x),L_1(x),L_2(x),L_3(x)\}$, by applying the Gram-Schmidt algorithm to the standard monomials $\{1,x,x^2,x^3\}$. The resulting Laguerre polynomials will form an orthogonal (or orthonormal if you include the normalization step in the Gram-Schmidt algorithm) family on $0<x<\infty$ with respect to the weight function $w(x)=e^{-x}$.
So, following the algorithm linked above (including the normalization) and using the weighted inner product above, you get (using the notation in the link):
\begin{align}
u_0(x)&=1,\\
u_1(x)&=x-{\langle x,1\rangle_w\over \langle 1,1\rangle_w}\cdot 1=x-{\int_0^\infty x\cdot 1\,e^{-x}\,dx\over \int_0^\infty 1^2\,e^{-x}\,dx}=x-1,\\
u_2(x)&=x^2-{\langle x^2,1\rangle_w\over \langle 1,1\rangle_w}\cdot 1-{\langle x^2,x-1\rangle_w\over \langle x-1,x-1\rangle_w}\cdot (x-1)\\
&=x^2-{\int_0^\infty x^2\cdot 1\,e^{-x}\,dx\over \int_0^\infty 1^2\,e^{-x}\,dx}\cdot 1-{\int_0^\infty x^2\cdot (x-1)\,e^{-x}\,dx\over \int_0^\infty (x-1)^2\,e^{-x}\,dx}\cdot (x-1)\\
&=x^2-4x+2,\\
u_3(x)&=x^3-{\langle x^3,1\rangle_w\over \langle 1,1\rangle_w}\cdot 1-{\langle x^3,x-1\rangle_w\over \langle x-1,x-1\rangle_w}\cdot (x-1)-{\langle x^3,u_2\rangle_w\over \langle u_2,u_2\rangle_w}\cdot (x^2-4x+2)\\
&=x^3-9x^2+18x-6.
\end{align}
To get the $L_k$, you just need to divide each $u_k$ above by its (weighted) norm:
\begin{align}
L_0(x)&={u_0\over \|u_0\|_w}={u_0\over \langle u_0,u_0\rangle_w^{1/2}}=1,\\
L_1(x)&={u_1\over \|u_1\|_w}={u_1\over \langle u_1,u_1\rangle_w^{1/2}}=x-1,\\
L_2(x)&={u_2\over \|u_2\|_w}={u_2\over \langle u_2,u_2\rangle_w^{1/2}}={x^2\over 2}-2x+1,\\
L_3(x)&={u_3\over \|u_3\|_w}={u_3\over \langle u_3,u_3\rangle_w^{1/2}}={x^3\over 6}-{3\over 2}x^2+3x-1.
\end{align}
and you get the first four Laguerre polynomials:
$$
\left\{1,x-1,\frac{x^2}{2}-2 x+1,\frac{x^3}{6}-\frac{3 x^2}{2}+3 x-1\right\}.
$$
Since these form an orthonormal family, you can expand $f(x)=e^{-2x}$ as follows:
$$
f(x)=\sum_{n=0}^k c_n L_n(x), \qquad c_n=\langle f,L_n\rangle_w=\int_0^\infty f(x)L_n(x)w(x)\,dx,
$$
where $k=0,1,2,3$ indicates (as you call it) the "degree" of the approximation, i.e., an expansion in terms of the first $k$ Laguerre polynomials.
For example,
\begin{align}
c_0&=\langle e^{-2x},1\rangle_w=\int_0^\infty e^{-2x}\cdot 1\cdot e^{-x}\,dx=1/3,\\
c_1&=\langle e^{-2x},x-1\rangle_w=\int_0^\infty e^{-2x}\cdot (x-1)\cdot e^{-x}\,dx=-2/9,\\
c_2&=\langle e^{-2x},\frac{x^2}{2}-2 x+1\rangle_w=\int_0^\infty e^{-2x}\cdot \left(\frac{x^2}{2}-2 x+1\right)\cdot e^{-x}\,dx=4/27,\\
c_3&=\langle e^{-2x},\frac{x^3}{6}-\frac{3 x^2}{2}+3 x-1\rangle_w=\int_0^\infty e^{-2x}\cdot \left(\frac{x^3}{6}-\frac{3 x^2}{2}+3 x-1\right)\cdot e^{-x}\,dx=-8/81.
\end{align}
Here are the $k=0,\dots,3$ approximations of $f(x)=e^{-2x}$:
Best Answer
Your calculations are correct.
The final result is
$P_3(x) = 0.0 + 0.5 (x-3) (x-2) (x+0)-4 (x-3) (x-1) (x+0)+4.5 (x-2) (x-1) (x+0) = x^3$
The formula for the error bound is given by:
$$E_n(x) = {f^{(n+1)}(\xi(x)) \over (n+1)!} \times (x-x_0)(x-x_1)...(x-x_n)$$
So, we have
$$E_3(x) = {f^{(4)}(\xi(x)) \over 4!} \times (x-0)(x-1)(x-2)(x-3)$$
The fourth derivative of $f(x) = x^3$ is zero, so $E_3(x) = 0$.
The reason for this is if $f(x) = $ polynomial of degree $M$ where $M \le N$, then $$f^{(n)}(x) = 0 \implies E_n(x) = 0 ~\forall~ x$$
Therefore $P_3(x)$ is an exact representation of $f(x) = x^3$.