The present study is related to the numerical solutions of new mathematical models based on the variable order Emden-Fowler delay differential equations. The shifted fractional Gegenbauer, $$C_{S,j}^{(\alpha ,\mu )}(t),$$
C
S
,
j
(
α
,
μ
)
(
t
)
,
operational matrices (OMs) of VO differentiation, in conjunction with the spectral collocation method are used to solve aforementioned models numerically. The VO operator of differentiation will be used in the Caputo sense. The proposed technique simplifies these models by reducing them to systems of algebraic equations that are easy to solve. To determine the effectiveness and accuracy of the sugested technique, the absolute errors and maximum absolute errors for four realistic models are studied and illustrated by several tables and graphs at different values of the VO and the SFG parameters; $$\alpha$$
α
and $$\mu.$$
μ
.
Also numerical comparisons between the suggested technique with other numerical methods in the existing literature are held. The numerical results confirm that the suggested technique is accurate, computationally efficient and easy to implement.