Here is a simple proof of the estimate that I gave in the comment. At the end I give some references.

We will use some potential-theoretic notions. Let $f(x) = g(x)$ in $B$ and $f(x) = h(x)$ in $C$. I understand that your notion of a "classical solution" asserts that
$$ u(x) = \int_{B \cup C} f(y) P_A(x, dy) , \tag{$\star$}$$
where $P_D$ stands for the harmonic measure for $(-\Delta)^s$ in an open set $D$. Note: formula ($\star$) holds true if, for example, $u$ is continuous and bounded on $A$ and $(-\Delta)^s u(x) = 0$ for every $x \in D$.

The expression for $P_D(x, dy)$ is known when $D$ is an interval, or a complement of the interval: we have
$$P_{(-r,r)}(x, dy) = c(s) \biggl(\frac{r^2 - x^2}{y^2 - r^2}\biggr)^s \frac{1}{|x - y|} \, dy,$$
and, provided that $s < \tfrac{1}{2}$,
$$P_{[-r,r]^c}(x, dy) = c(s) \biggl(\frac{x^2 - r^2}{r^2 - y^2}\biggr)^s \frac{1}{|x - y|} \, dy$$
(a similar expression is available for $s \geqslant \tfrac{1}{2}$ as well, but it is slightly more complicated). The first one is given (in an arbitrary dimension) in Riesz's 1938 paper, and I think the other one appeared first in the paper of Blumenthal, Getoor and Ray from 1961. (In dimension 1, however, these expressions may be even older, I think Weyl considered similar problems; however, I do not have a reference at hand).

It is well-known that $P_D(x, dy)$ is domain-monotone: $P_{D_1}(x, E) \leqslant P_{D_2}(x, E)$ if $D_1 \subseteq D_2$ and $D_2 \cap E = 0$. Thus,
$$ P_A(x, dy) \leqslant P_{(-2,2)}(x, dy) \qquad \text{on $C$,} $$
and
$$ P_A(x, dy) \leqslant P_{[-1,1]^c}(x, dy) \qquad \text{on $B$.} $$
By combining these estimates with explicit expressions listed above and plugging them into ($\star$), we get
$$ |u(x)| \leqslant c(s) \int_C \biggl(\frac{4 - x^2}{y^2 - 4}\biggr)^s \frac{1}{|x - y|} \, |h(y)| dy + c(s) \int_B \biggl(\frac{x^2 - 1}{1 - y^2}\biggr)^s \frac{1}{|x - y|} \, |g(y)| dy.$$

Some references: all the information required to find the above estimate are given already in:

M. Riesz, Intégrales de Riemann–Liouville et potentiels, Acta Sci. Math. Szeged, 9 (1938), 1–42

and in:

R. M. Blumenthal, R. K. Getoor, and D. B. Ray, On the distribution of first hits for the symmetric stable processes, Trans. Am. Math. Soc., 99 (1961), 540–554

Similar bounds in higher dimensions are related to the boundary Harnack inequality for $(-\Delta)^s$, a business started in

K. Bogdan, The boundary Harnack principle for the fractional Laplacian, Stud. Math., 123(1) (1997), 43–80.

Most papers, however, consider dimension at least $2$. I believe there is a reference that discusses dimension $1$ explicitly, but I cannot recall it right now (actually, my paper with Bogdan and Kulczycki has it hidden somewhere, I think). The case $s \geqslant \tfrac{1}{2}$ is more problematic for the above approach, but not really for the boundary Harnack inequality methods. In the end one gets a bound of the form
$$ |u(x)| \leqslant c(A, s) \int_{A^c} \biggl(\frac{\delta_A(x)}{\delta_A(y)}\biggr)^s \frac{1}{|x - y|} \, |f(y)| dy ,$$
no matter what the value of $s \in (0, 1)$ is, where $\delta_A(x)$ is the distance of $x$ from the boundary of $A$. Here $A$ can be an arbitrary bounded $C^{1,1}$ set in $\mathbb{R}^n$.

For more information, see my survey:

M. Kwaśnicki, *Fractional Laplace Operator and its Properties*, in: A. Kochubei, Y. Luchko, *Handbook of Fractional Calculus with Applications. Volume 1: Basic Theory*, De Gruyter Reference, De Gruyter, Berlin, 2019