  We consider the homogeneous equation ${\mathcal A} u=0$, where ${\mathcal A}$
is a symmetric and coercive elliptic operator in $H^1(\Omega)$ with $\Omega$
bounded domain in ${{\mathbb R}}^d$. The boundary conditions involve fractional
power $\alpha$, $ 0 < \alpha <1$, of the Steklov spectral operator arising in
Dirichlet to Neumann map. For such problems we discuss two different numerical
methods: (1) a computational algorithm based on an approximation of the
integral representation of the fractional power of the operator and (2)
numerical technique involving an auxiliary Cauchy problem for an
ultra-parabolic equation and its subsequent approximation by a time stepping
technique. For both methods we present numerical experiment for a model
two-dimensional problem that demonstrate the accuracy, efficiency, and
stability of the algorithms.
