I estimated (w.r.t. $\varepsilon$) the expression
\begin{align} &\left|\int_{-1}^{x_0-\varepsilon} (1-x)^{n-p}(1+x)^p+\int_{x_0+\varepsilon}^1 (1-x)^{n-p}(1+x)^p \, dx \right | \\[6pt] \leqslant {} & \left| (x_0-\varepsilon+1)(1-x_0+\varepsilon)^{n-p}(1+x_0-\varepsilon)^p+(1-x_0-\varepsilon)(1-x_0-\varepsilon)^{n-p}(1+x_0+\varepsilon)^{p} \right| \\[6pt] = {} & \left|(1-x_0+\varepsilon)^{n-p}(1+x_0-\varepsilon)^{p+1}+(1-x_0-\varepsilon)^{n-p+1}(1+x_0+\varepsilon)^p \right| \end{align}
where $x_0={2p-n\over n}$ is the maximum of the integrand $(1-x)^{n-p}(1+x)^{p}$ in the interval $[-1,1]$ and $p\le n$ and $n,p \in \Bbb N $.
This is possible since the function goes to zero monotonously from $x_0$.
Does anyone know a way to make a better estimate?