In this paper, a vaccination model for SARS-CoV-2 variants is proposed and is studied using fractional differential operators involving a non-singular kernel. It is worth mentioning that variability in transmission rates occurs because of the particular population that is vaccinated, and hence, the asymptomatic infected classes are classified on the basis of their vaccination history. Using the Banach contraction principle and the Arzela–Ascoli theorem, existence and uniqueness results for the proposed model are presented. Two different numerical approaches, the fractional Euler and Lagrange polynomial methods, are employed to approximate the model’s solution. The model is then fitted to data associated with COVID-19 deaths in Pakistan between 1 January 2022 and 10 April 2022. It is concluded that our model is much aligned with the data when the order of the fractional derivative ζ=0.96. The two different approaches are then compared with different step sizes. It is observed that they behave alike for small step sizes and exhibit different behaviour for larger step sizes. Based on the numerical assessment of the model presented herein, the impact of vaccination and the fractional order are highlighted. It is also noted that vaccination could remarkably decrease the spikes of different emerging variants of SARS-CoV-2 within the population.