We consider the homogeneous equation Au = 0, where A is a symmetric and coercive elliptic operator in H1(Ω) with Ω bounded domain in ℝd. The boundary conditions involve fractional power α, 0 < α < 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. © 2017 Diogenes Co., Sofia.