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

 A: Here is a solution:
$$
\frac{1}{2}+\frac{1}{\pi}(\frac{x}{\epsilon}-\frac{x^3}{3\epsilon^{3}}+\frac{x^{5}}{5\epsilon^5}\cdots)
$$
A: It is not possible to use polynomial as Heaviside step function with a good average precision, because any polynomial is infinite at both positive and negative infinity and Heaviside is not. Second, you would need to have zero value for all negative values. Even if you take that the value is only approximately zero you would need all larger and larger degree to achieve that.
In that sense the only option is to use rational function of two polynomials.
If you still insist of having just one polynomial expression then you can use normal Lagrange/Newton approximations taking as many points as you want. This can give you some accuracy within some region, and for the jump around $0$ you need to take more points around $0$.
Here is the polynomial for:
$$(0,1/2) (-1,0) (-2,0) (-3,0) (-4,0) (1,1) (2,1) (3,1) (4,1) (5,1)$$
$$\frac{x^9}{51840}-\frac{13x^7}{12096}+\frac{71x^5}{3456}-\frac{107x^3}{648}+\frac{1627x}{2520}+\frac{1}{2}$$
Is it possible to get the coefficients directly?
Sure, just create a linear system and solve it against the unknown coefficients. Notice that there is an infinite polynomial given by Taylor expansion of any of the approximations like the one using $\tanh(kx)$ with extremely slow convergence.
In essence, since we ask a polynomial to be too much flat, a polynomial is a bad approximation to step function no matter what we do. A rational polynomial function is much better choice.
A: 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)).$$
