  We consider a least-squares variational kernel-based method for numerical
solution of second order elliptic partial differential equations on a
multi-dimensional domain. In this setting it is not assumed that the
differential operator is self-adjoint or positive definite as it should be in
the Rayleigh-Ritz setting. However, the new scheme leads to a symmetric and
positive definite algebraic system of equations. Moreover, the resulting method
does not rely on certain subspaces satisfying the boundary conditions. The
trial space for discretization is provided via standard kernels that reproduce
the Sobolev spaces as their native spaces. The error analysis of the method is
given, but it is partly subjected to an inverse inequality on the boundary
which is still an open problem. The condition number of the final linear system
is approximated in terms of the smoothness of the kernel and the discretization
quality. Finally, the results of some computational experiments support the
theoretical error bounds.
