The approximate analytical expressions of tripole-mode and quadrupole-mode solitons in (1 + 1)-dimensional nematic liquid crystals are obtained by applying the variational approach. It is found that the soliton powers for the two types of solitons are not equal with the same parameters, which is much different from their counterparts in the Snyder-Mitchell model (an ideal and typical strongly nolocal nonlinear model). The numerical simulations show that for the strongly nonlocal case, by expanding the response function to the second order, the approximate soliton solutions are in good agreement with the numerical results. Furthermore, by expanding the respond function to the higher orders, the accuracy and the validity range of the approximate soliton solutions increase. If the response function is expanded to the tenth order, the approximate solutions are still valid for the general nonlocal case.