In this paper, a nonlinear fractional order model of COVID-19 is approximated. For this aim, at first we apply the Caputo–Fabrizio fractional derivative to model the usual form of the phenomenon. In order to show the existence of a solution, the Banach fixed point theorem and the Picard–Lindelof approach are used. Additionally, the stability analysis is discussed using the fixed point theorem. The model is approximated based on Indian data and using the homotopy analysis transform method (HATM), which is among the most famous, flexible and applicable semi-analytical methods. After that, the CESTAC (Controle et Estimation Stochastique des Arrondis de Calculs) method and the CADNA (Control of Accuracy and Debugging for Numerical Applications) library, which are based on discrete stochastic arithmetic (DSA), are applied to validate the numerical results of the HATM. Additionally, the stopping condition in the numerical algorithm is based on two successive approximations and the main theorem of the CESTAC method can aid us analytically to apply the new terminations criterion instead of the usual absolute error that we use in the floating-point arithmetic (FPA). Finding the optimal approximations and the optimal iteration of the HATM to solve the nonlinear fractional order model of COVID-19 are the main novelties of this study.