**1. Introduction**

Systems of ordinary differential equations are useful in representing some real life problems in terms of the mathematical expressions, which abound in the fields of biological, physical, engineering, financial or sociological fields. It is well known that many nonlinear problems in these fields can be well modeled by systems of ordinary differential equations. However, finding exact solutions of systems of ordinary differential equations involving nonlinear terms can be extremely difficult in most of the situations. In addition, we know that exact solutions of most realistic systems of ordinary differential equations cannot be found, so we need numerical and approximate methods for finding approximate solutions.There are a lot of methods that have been studied by many researchers to solve the systems of ordinary differential equations. Some of these methods are the multi-step method proposed by Hojjati et al. [1], the collocation method presented by Mastorakis [2], the Adomian decomposition method improves [3], the exponential Galerkin method introduced by Yüzba¸sı and Karaçayır [4], the exponential collocation method proposed by Yüzba¸sı [5], the Galerkin finite element method given by Al-Omari et al. [6].

In this study, we are interested in the SIR model, a model of an epidemic of an infectious disease in a population. This model comprises three types of individuals: those who might be susceptible to the disease, those who might be infected with the disease , and those who might have recovered or be immune from the disease. The model thus has three classes or states.

The following system determines the progress of the disease [7]:

$$\begin{aligned} \frac{dS}{dt} &= -\beta S(t)I(t) \\ \frac{dI}{dt} &= \beta S(t)I(t) - \gamma I(t) \\ \frac{dR}{dt} &= \gamma I(t) \end{aligned} \tag{1}$$

with initial conditions

$$S(0) = N\_{\mathbb{S}\prime}I(0) = N\_{I\prime}R(0) = N\_{\mathbb{R}}.\tag{2}$$

*NS* = the number of susceptible individuals in the population at time *t*.

*NI* = the number of infected individuals in the population at time *t*.

*NR* = the number of recovered individuals in the population at time *t*.

*N* = the population size.

*β* = the transmissivity rate.

*γ* = is the recovery rate. Note that, at any given time, an individual can only be in one of the three groups. Thus, *NS* + *NI* + *NR* = *N*.

Finding exact solutions of SIR models is important because biologists could use it to design and run experiments to observe the spread of infectious diseases by introducing natural initial conditions. Through these experiments, as well as through mathematical modelling, one can learn the ways on how to control the spread of epidemics. It is extremely difficulty to obtain the exact solutions for such problems that actually represented such phenomena. It is a big task for scientific community to search for appropriate methods. Within two decades, to obtain approximate solutions of Equation (1), some authors have studied this model using different methods. For Equation (1), Argub and El-Ajou used Homotopy Analysis Method for different parameter values [7], Awawdeh et al. used Homotopy Analysis Method [8], Biazar used the Adomian decomposition method [9], Rafei et al. applied homotopy perturbation method [10], Ibrahim et al. applied Differential Transformation Approach [11]. In [12], this system was solved using Laplace-Adomian decomposition method. In [13], Harman and Johnston solved the epidemic model using stochastic Galerkin method. Equation (1) was solved using 4th order Runge-Kutta method by Kousar et al. [14] and using Euler, Runge Kutta-2 and Runge-Kutta-4 methods by Hussain et al. [15].

The collocation method has become progressively favourite to solve differential equations. This method can reduce the complexity of solving the systems of ordinary differential equations for epidemic models with high dimensions and it is very useful in contributing highly accurate solutions to differential equations. In this study, Hermite polynomials, a class of the orthogonal polynomials {*H*0(*t*), *H*1(*t*), ... , *HL*(*t*)} that are orthogonol on (−∞, ∞), are used. Hermite polynomials have advantages over other orthogonal polynomials. Hermite collocation method (HCM) has been used to solve systems of nonlinear ordinary differential equations with special initial conditions. The most important advantage of the presented method is that it transforms this system (1) into a nonlinear system of algebraic equations which can be easily solved. Until recently, HCM has been used to obtain solutions to a higher-order linear Fredholm integro differential equations in [16], to linear fractional order Systems of differential equations in [17], to differential difference equations in [18], to fractional order differential equations in [19] and to the neutral functional-differential equations with proportional delays in [20].
