Optical spatial solitary waves are considered in a nonlocal thermal focusing medium with non-symmetric boundary conditions. The governing equations consist of a nonlinear Schrödinger equation for the light beam and a Poisson equation for the temperature of the medium. Three numerical methods are investigated for calculating the ground and excited solitary wave solutions of the coupled system. It is found that the Newton conjugate gradient method is the most computationally efficient and versatile numerical technique. The solutions show that by varying the ambient temperature, the solitary wave is deflected towards the warmer boundary. Solitary wave stability is also examined both theoretically and numerically, via power versus propagation constant curves and numerical simulations of the governing partial differential equations. Both the ground and excited state solitary waves are found to be stable. The Newton conjugate gradient method should also prove extremely useful for calculating solitary waves of other related optical systems, which support nonlocal spatial solitary waves, such as nematic liquid crystals. © 2013 IOP Publishing Ltd.