Could any one tell me how to solve this one?
Given $f\in C[0,\infty)$ such that $f(x)\to 0$ as $x\to\infty$ we need to show that for any $\epsilon>0$ there is a polynomial $p$ such that $|f(x)-e^{-x}p(x)|<\epsilon \qquad \forall~ x\in[0,\infty)$
I just know the statement of Weierstrass Polynomial Approximation Theorem and that seems very far from the given problem, but somehow I feel I need to apply the theorem.
It's not a "solution" (proof?) but I would like to give some details to your question.
$$ \vert f(x)-e^{-x}p(x)\vert=e^{-x}\vert g(x)-p(x)\vert, $$ where $\lim_{x\to\infty}e^{-x}g(x)=0$. Your question is about the density of polynomials in the weighted space (weight function is $e^{-x}$) on the half-line. More info can be found Stone-Weierstrass theorem (locally compact version). The literature of weighted approximation in $L_p$ spaces is huge.
We may consider $f$ defined on the compact space $[0, \infty]$, continuous, $f(\infty) =0$.
Consider the function $\phi(x) = \exp(-x)$. It separates the points and its only $0$ point is $\infty$. From Stone-Weierstrass theorem we conclude that the uniform closure $\mathbb{R}$- algebra of functions $\sum_{k\ge 1} a_k e^{-k x}$ is the algebra of continuous functions on $[0, \infty]$ that are $0$ at $\infty$.
Now, we want to approximate every $e^{-(k+1) x}$ ($k\ge 0$) by some $e^{-x}P(x)$, with $P(x)$ a polynomial.
Let's try first to approximate $e^{-2x}$. To do this we approximate $e^{-x}$ with its Taylor series and multiply. Note that we have $$\sum_{k=0}^{2m} (-1)^k \frac{x^k}{k!} > e^{-x} > \sum_{k=0}^{2n+1} (-1)^k \frac{x^k}{k!}$$ for $x>0$ for any $m$, $n$ natural numbers. We conclude that $$|e^{-x} - T_{n-1}(x)| < \frac{x^n}{n!}$$ where $T_{n-1}$ is the Taylor polynomial of order $n-1$ for $e^{-x}$. From here we get $$|e^{-2x} - T_{n-1} e^{-x}| < e^{-x} \frac{x^n}{n!}$$
Now it is left to show that $\max_{x\ge 0} e^{-x} \frac{x^n}{n!} \to 0$ as $n\to \infty$. It is easy to see that for each $n$ the corresponding $\max$ is at $x=n$, so we need to see that $$\frac{\left(\frac{n}{e}\right)^n}{n!}\to 0$$ as $n\to \infty$. Now Stirling formula is your friend.
One would be tempted to do similarly for $e^{-3x}$, $e^{-4x}$ and be done with it. However, the Taylor approximation estimate breaks down rather badly: the corresponding maximum will diverge as $n\to \infty$.
Another approach: try to show that the closure of $(e^{-x} P(x))_P$ is closed under multiplication. That would do it. So it's enough to show that $e^{-x} P(x) e^{-x} Q(x)$ can be approximated as $e^{-x} R(x)$. So it would be enough to show that $e^{-2x} x^l$ can be approximated for all $l\ge 0$. Again, Taylor expansion does not seem to work.
In the end, we recall the Laguerre polynomials, and the fact that the associated functions form an orthogonal basis in $L^2[0, \infty)$. So now at least we know how to get for a nice function an approximating $e^{-x} P(x)$. This actually seems to be working very well. An equivalent approach is to consider Hermite polynomials on $\mathbb{R}$. The advantage is that the formulas seem to be neater. We can get ( at least formally) the expansion of $e^{-x^2/(1-\rho^2)}$ in terms of $H_{2n}(x) e^{-x^2}$ by using the Mehler formula $$\frac{\exp( - \frac{x^2 - 2 \rho x y + y^2}{1-\rho^2})}{\sqrt{1-\rho^2}} = \sum_{n\ge 0} \frac{(u/2)^n}{n!} H_n(x) e^{-x^2} H_n(y) e^{-y^2}$$
( see a WA testing). Plug in the above $y=0$ and we get a series expansion $$e^{-x^2/(1-\rho^2)} = \sqrt{1-\rho^2} \sum_{n\ge 0} \frac{(\rho/2)^n}{n!}H_n(0) H_n(x) e^{-x^2}$$
(see with WA). Note that $H_n(0) = 0$ for $n$ odd and $H_{2n}(x)$ is an even polynomial. Also, these are the physicists' polynomials.
We can wave our hands and claim that we have uniform convergence over $\mathbb{R}$ in the above expansion. This is probably so, but I don't have a source for that. More details are needed.
$\bf{Added:}$ It turns out that we can use Laguerre polynomials and the arguments are simpler. Consider the expansion
$$\frac{\exp(-\frac{x t}{1-t})}{1-t} = \sum_{n=0}^{\infty} L_n(x) t^n$$ where we can consider this just a formal expansion in $t$. Multiplying on both sides by $\exp(-x/2)$ we get $$\exp(-\frac{1+t}{2(1-t)} x ) = (1-t)\sum_{n\ge 0} t^n L_n(x) e^{-x/2} $$
We look at it now as an expansion of a function of $x$, with $t$ a parameter. We can rewrite it as
$$e^{-a x} = \frac{2}{2a+1}\sum_{n\ge 0} (\frac{2a-1}{2a+1})^n L_n(x) e^{-x/2}$$
Now, if $\operatorname{Re} a>0$ then $|\frac{2a-1}{2a+1}|<1$. Important fact: $\sup_{x\ge 0} |L_n(x)| e^{-x/2} = L_n(0) = 1$ for all $n$. ( reference here ). We conclude that in the above we have normal convergence in $x$. ( $a$ in a compact subset of the right half plane). In particular every $e^{-a x}$ can be approximated uniformly by functions $e^{-x/2} P(x)$.
The functions $L_n(x) e^{-x/2}$ form an orthonormal system in $L^2[0, \infty)$. The Laguerre expansion of every function $e^{-a x}$ converges normally to $e^{-a x}$. We could obtain the expansion of other functions if we knew their inverse Laplace transform.
