How to approximate Heaviside function by polynomial I have a Heaviside smooth function that defined as
$$H_{\epsilon}=\frac {1}{2} [1+\frac {2}{\pi} \arctan(\frac {x}{\epsilon})]$$

I want to use polynominal to approximate the Heaviside function. Could you suggest to me a solution? Thanks
UPDATE: This is Bombyx mori result in blue line and my expected result is red line


Althogh I think Bombyx mori's solution is great and simple, maybe the following approximation is better for you on the interval $[-1,1]$.
$$
p_\varepsilon(x) = -\tfrac{1}{6\pi}\left(16\arctan\left(\tfrac{1}{2\varepsilon}\right)-8\arctan\left(\tfrac{1}{\varepsilon}\right)\right)x^3 - \tfrac{1}{6\pi}\left(2\arctan\left(\tfrac{1}{\varepsilon}\right)-16\arctan\left(\tfrac{1}{2\varepsilon}\right)\right)x+\tfrac12.
$$
I've got this by making a polynomial interpolation for the data points: $$(-1,H_\varepsilon(-1)), (-1/2,H_\varepsilon(-1/2)), (0,H_\varepsilon(0)),(1/2,H_\varepsilon(1/2)), (1,H_\varepsilon(1)).$$
