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.