  The nonlinear Schr\"{o}dinger equation with a linear periodic potential and a
nonlinearity coefficient $\Gamma$ with a discontinuity supports stationary
localized solitary waves with frequencies inside spectral gaps, so called
surface gap solitons (SGSs). We compute families of 1D SGSs using the arclength
continuation method for a range of values of the jump in $\Gamma$. Using
asymptotics, we show that when the frequency parameter converges to the
bifurcation gap edge, the size of the allowed jump in $\Gamma$ converges to 0
for SGSs centered at any $x_c\in \R$.
  Linear stability of SGSs is next determined via the numerical Evans function
method, in which the stable and unstable manifolds corresponding to the 0
solution of the linearized spectral ODE problem need to be evolved. Zeros of
the Evans function coincide with eigenvalues of the linearized operator. Far
from the SGS center the manifolds are spanned by exponentially
decaying/increasing Bloch functions. Evolution of the manifolds suffers from
stiffness but a numerically stable formulation is possible in the exterior
algebra formulation and with the use of Grassmanian preserving ODE integrators.
Eigenvalues with positive real part above a small constant are then detected
using the complex argument principle and a contour parallel to the imaginary
axis. The location of real eigenvalues is found via a straightforward
evaluation of the Evans function along the real axis and several complex
eigenvalues are located using M\"{u}ller's method. The numerical Evans function
method is described in detail. Our results show the existence of both unstable
and stable SGSs (possibly with a weak instability), where stability is obtained
even for some SGSs centered in the domain half with the less focusing
nonlinearity. Direct simulations of the PDE for selected SGS examples confirm
the results of Evans function computations.
