Approximating step functions with polynomials Let $t_1 < t_2 < \cdots <t_m$ be real, and $X = \cup_{i=1}^{m-1} (t_i, t_{i+1})$ be a union of real open intervals.  Let $f:X \rightarrow \{-1, 1\}$ be any piecewise constant function of form
$$
f(x) = 
    \begin{cases} 
      a_1 & \text{ if } t_1 < x < t_2 \\
      a_2 & \text{ if }t_2 < x < t_3 \\
      \vdots \\
a_{m-2} & \text{ if } t_{m-2} < x < t_{m-1} \\ 
      a_{m-1} & \text{ if } t_{m-1} < x < t_m 
   \end{cases}
$$
Where $a_i \in \{-1, 1\}$, and $a_{i} = -a_{i+1}$ for $i = 1, ..., m-1$.
I have a number of questions regarding polynomial approximations of such a function $f$:

*

*Can we always find a sequence of polynomials $(p_n)$ so that $(p_n)$ converge pointwise to $f$, AND we have some fixed (not arbitrary) global error bound, say $1$, such that $|p_n(x) - f(x)| \leq 1$ for all $x \in X$ and $n \in \mathbb{N}$?

*If so, are such polynomials easy to find?

*How quickly do we get convergence?

I am aware that, upon picking a suitable inner product, we can use any collection of orthonormal polynomials to make approximations of functions.  For example I know the Chebyshev, Bernstein, Jacobi etc. polynomials can be used to approximate continuous functions on bounded intervals, but I have found no theorem that says we can use these to construct approximations for arbitrary piecewise constant functions like the one given above.
Indeed, it is easy to find a polynomial approximation for the Heaviside Step function for example, however it is unclear how, or if this an be done for more complicated step functions.

The given step function can be presented in the form of
$$f(x)=(-1)^m a^\ _1\operatorname{sgn}y(x),\quad\text{where}\quad y(x)=\prod_{j=1}^{m-1}(x-t_j).\tag1$$
Function $\,\operatorname{sgn}z\,$ can be presented via Fourier series in the form of
$$\operatorname{sgn}z=\dfrac 4\pi\sum_{n=1}^\infty \dfrac{\sin((2n-1)\pi z)}{2n-1},\tag2$$

(n=1...20)
where the sine functions are the polynomials of $\,\sin \pi z:\,$
$$\sin((2n-1)\pi z)=P_{2n-1}(\sin \pi z),\tag3$$
$$P_{2n-1}(t) = t\sum_{j=0}^{n-1} (-1)^j\dbinom{2n-2-j}j(4-4t^2)^{n-1-j}.\tag4$$
Therefore,
$$\operatorname{sgn}z=\dfrac4\pi\,\lim\limits_{N\to\infty} \sum_{n=1}^N \dfrac1{2n-1}P_{2n-1}(\sin \pi z).\tag5$$
Since the function $\operatorname{sgn} z$ has the gap in the point $\,z=0,\,$ then attempt of transformation of formulas $(5)$ to the polynomial series leads to the divergent expressions for the coefficients.
However, looks possible to build the polynomial approximation of the function $f(x)$ as the polynomial sequence $\,f_N(x),\,$ based on the formulas $(1),(4),(5)$ and the known Maclaurin series for $\,\sin \pi z.\,$
The optimization task solution (via logarithmic derivative) gives
$$Y=\max\limits_{x\in(t_1,t_m)} |y(x)|=\max(|y(x_l)|,|y(x_h)|),\tag6$$
where
$$\dfrac1{x_l-t_1}=\prod_{j=2}^{m-1}\dfrac1{t_j-x_l},\quad x_l\in(t_1,t_2)\tag7$$
and
$$\dfrac1{x_h-t_{m-1}} = \prod_{j=1}^{m-2}\dfrac1{x_h-t_j},\quad x_h\in(t_{m-1}, t_m).\tag8$$
Then
$$\operatorname{sgn}(y(x)) = \operatorname{sgn} z(x),\quad\text{where}\quad z(x)=\dfrac {y(x)}Y,\quad |z(x)|\le1.\tag9$$
Applying of $\,z(x)\,$ instead $\,y(x)\,$ in $(1)$ allows to control the calculation accuracy of the sine function and finally - the required distances of the given function gaps, which provide the required calculation accuracy of the approximations $f_N(x)$.
