In analogy to a perturbed harmonic oscillator, we calculate the fundamental and some other higher order soliton solutions of the nonlocal nonlinear Schödinger equation (NNLSE) in the second approximation in the generally nonlocal case. Comparing with numerical simulations we show that soliton solutions in the second approximation can describe the generally nonlocal soliton states of the NNLSE more exactly than that in the zeroth approximation. We show that for the nonlocal case of an exponential-decay type nonlocal response the Gaussian-function-like soliton solutions cannot describe the nonlocal soliton states exactly even in the strongly nonlocal case. The properties of such nonlocal solitons are investigated. In the strongly nonlocal limit, the soliton’s power and phase constant are both in inverse proportion to the fourth power of its beam width for the nonlocal case of a Gaussian function type nonlocal response, and are both in inverse proportion to the third power of its beam width for the nonlocal case of an exponential-decay type nonlocal response.