Here is a solution in an integral form. I suppose it can be written in terms of hypergeometric functions (or maybe Meijer G-functions), but I did not attempt to do that. Once this is done, extension to general $\Re k > -1$ should follow by analytic continuation.
For $k \in \mathbb{C}$ let
$$ f_k(x) = \begin{cases} x^k & \text{for $x > 0$,} \\ 0 & \text{otherwise.} \end{cases} $$
Let $L = (-\Delta)^{s/2} = |\nabla|^s$ denote the fractional Laplacian.
Lemma: If $-1 < \Re k < s$, we have
$$ L f_k(x) = \begin{cases} a_k x^{k - s} & \text{for $x > 0$,} \\ b_k (-x)^{k - s} & \text{for $x < 0$,} \end{cases} $$
where
$$ a_k = 2^{s-1} \frac{\Gamma(\tfrac{k+1}{2})\Gamma(\tfrac{k}{2})}{\Gamma(\tfrac{k+1-s}{2})\Gamma(-\tfrac{k}{2})} + 2^{s-1} \frac{\Gamma(\tfrac{k}{2}+1)\Gamma(\tfrac{k-1}{2})}{\Gamma(\tfrac{k-s}{2}+1)\Gamma(-\tfrac{k-1}{2})} $$
and
$$ a_k = 2^{s-1} \frac{\Gamma(\tfrac{k+1}{2})\Gamma(\tfrac{k}{2})}{\Gamma(\tfrac{k+1-s}{2})\Gamma(-\tfrac{k}{2})} - 2^{s-1} \frac{\Gamma(\tfrac{k}{2}+1)\Gamma(\tfrac{k-1}{2})}{\Gamma(\tfrac{k-s}{2}+1)\Gamma(-\tfrac{k-1}{2})} \, . $$
Proof: In $\mathbb{R}^n$, it is known that
$$ L[|x|^k] = 2^s \frac{\Gamma(\tfrac{k+n}{2})\Gamma(\tfrac{k}{2})}{\Gamma(\tfrac{k+n-s}{2})\Gamma(-\tfrac{k}{2})} \, |x|^{k - s} $$
and
$$ L[|x|^{k - 1} x_1] = 2^s \frac{\Gamma(\tfrac{(k-1)+(n+2)}{2})\Gamma(\tfrac{k-1}{2})}{\Gamma(\tfrac{(k-1)+(n+2)-s}{2})\Gamma(-\tfrac{k-1}{2})} \, |x|^{(k - 1) - s} x_1 ; $$
the first identity is quite standard, the latter one is also likely well-known, and both follow, for example, from Theorem 1 in my paper Fractional Laplace operator and Meijer G-function with Bartłomiej Dyda and Alexey Kuznetsov, or Theorem 3.6 in my survey Fractional Laplace Operator and its Properties. Taking $n = 1$ and combining both identities, we get the desired result. $\square$
Corollary: Let $-1 < \Re k < s$,
$$ v_k(x) = \frac{1}{a_{k+s}} \, f_{k+s}(x) - \frac{1}{a_{k+s} \Gamma(1 + \tfrac{s}{2}) |\Gamma(-\tfrac{s}{2})|} \int_1^\infty \frac{(x - x^2)^{s/2}}{(y^2 - y)^{s/2} (y - x)} \, f_{k+s}(y) dy $$
for $x \in (0, 1)$, and $v_k(x) = 0$ otherwise. Then $L v_k(x) = f_k(x)$ for $x \in (0, 1)$.
Proof: Note that $v_k$ is a difference of $f_{k+s} / a_{k+s}$ and an $L$-harmonic function in $(0, 1)$ (the integral term in the definition is just the $L$-harmonic reduction of $f_{k+s} / a_{k+s}$, that is, the integral of $f_{k+s} / a_{k+s}$ with respect to the Poisson kernel for $L$). Thus, $L v_k = L(f_{k+s}) / a_{k+s}) = f_k$ in $(0, 1)$ by our lemma. $\square$
The above works for any $s$ such that $\Re s > -1$, I suppose.
