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}$:
Let's call the first 3 orthogonal polynomials we are looking for $v_0(x), v_1(x), v_2(x)$. We are going to solve for $v_0(x)$ by starting with $u_0(x) = 1$ as your first polynomial. You can find the normalization constant $N$ by solving:
$\int^\infty_0 v_0(x)^2 e^{-x}dx = \frac{1}{N^2} $
Your normalized polynomial $v_0(x)$ will be $Nu_0(x)$. In this case, you will find that N=1, or that your function $v_0(x) = 1$ is already normalized.
Then continue with $u_1(x) = x$. You have to make it orthogonal with your first orthogonal polynomial $v_0(x)$ by evaluating the integral:
$\int^\infty_0 u_1(x)v_0(x) e^{-x}dx$
You have to subtract the result of this integral from $u_1(x)$ to get your first orthogonal polynomial. The integral turns out to be -1 (solved by integration by parts), and so you get $v_1 = x - 1$. You can find that the normalization constant is 1 by doing the integration
$\int^\infty_0 v_1(x)^2 e^{-x}dx = \frac{1}{N^2} $
Finally,to find $v_2(x)$, we start with $u_2(x) = x^2$ and make it orthogonal to both of our polynomials $v_0(x)$ and $v_1(x)$. This time we have to subtract the following two integrals from $u_2(x)$:
$\int^\infty_0 u_2(x)v_0(x) e^{-x}dx$
$\int^\infty_0 u_2(x)v_1(x) e^{-x}dx$.
This should lead to $v_2(x) = x^2 - 4x + 2$. Do the usual normalization integral, and this time you will find that $N=\frac{1}{2}$.
You can check your answers here: https://en.wikipedia.org/wiki/Laguerre_polynomials#The_first_few_polynomials
(Note that our $v_1(x)$ has a minus sign as compared to theirs. That is fine, since a minus sign doesn't change the fact the orthnormality of the polynomial)
Best Answer
We apply the Gram Schmidt orthogonalization procedure to $$ P_1(x)=c \qquad P_2(x)=ax \qquad P_3(x)=bx^2$$
1) Take $P_1$ and normalize it
$$ ||P_1(x)||^2=\int_0^{\infty}c^2e^{-x}=c^2 $$
So $c=1 $ and $L_0=1 $
2) Take $P_2$ and subtract to it its projection on $L_0$, then normalize the result
$$P_{2|L_0} = \int_0^{\infty} (ax)e^{-x}dx =a \\P_2-P_{2|L_0}=a(x-1) \\ ||P_2-P_{2|L_0}||^2 = \int_0^{\infty} a^2(x-1)^2e^{-x}=a^2$$
Now we are free to choose either $ a=1$ or $a=-1 $ , to remain consistent with your notation we choose the latter and get
$$L_1=1-x$$
3)Take $P_3$ and subtract to it its projection on $L_0$ and on $L_1$ then normalize it and get $L_2$ with very analogous procedure