I have a strong impression that something like that has been asked before (perhaps, by somebody else) but it is easier to answer again than to find that old thread.
What you ask for is patently impossible. Indeed, assume that you have a polynomial $f_n(z)$ that approximates $\sqrt z$ on the interval $[\frac 13,1]$ with precision $e^{-cn}$. Consider the domain $\Omega=\{z: \frac 13\le |z|\le 1, 0\le \arg z\le \frac{3\pi}2\}$.
The function $g(z)=f_n(z)-\sqrt{z}$ is then analytic in $\Omega$, continuous up to the boundary, and bounded by $B_n+1$. By the standard two constant lemma, we have $$ |g(-2/3)|\le [\max_{\partial\Omega\setminus[\frac 13,1]}|g|]^{1-\gamma}[\max_{[\frac 13,1]}|g|]^{\gamma}\le (B_n+1)^{1-\gamma}e^{-\gamma cn} $$ with some constant $\gamma\in(0,1)$ (the harmonic measure of $[\frac 13,1]$ with respect to the domain $\Omega$ and the point $-2/3$). If $B_n$ is subexponential in $n$, the RHS tends to $0$ as $n\to\infty$, so we get $f_n(-2/3)$ close to $i\sqrt{\frac 23}$.
Considering the domain symmetric to $\Omega$ with respect to the real axis, we conclude that $f_n(-2/3)$ must be also very close to $-i\sqrt{\frac 23}$. But those two numbers are rather far apart.
