The reduced Lagrange density of nonlocal nonlinear Schrdinger equation(NNLSE) is obtained by expanding the real symmetric response function in Taylor’s series in strongly nonlocal Kerr media. The problem of higher-order beam propagation can be analyzed by a variational approach, the equations are obtained for the evolution during propagation of the parameters of the trial solution and exact analytical Hermite-Gaussian(HG) solutions are found. HG solitons are formed when the input power is equal to the critical power. We demonstrated that the analytical HG solutions are in good agreement with the numerical simulations in the case of strong nonlocality.