The purpose of this paper is to study the existence and multiplicity of positive solutions for a mathematical model of thermal explosion which is described by the system$$\left\{\begin{array}{ll}-\Delta u = \lambda f(v), & x\in \Omega,\\-\Delta v = \lambda g(u), & x\in \Omega,\\\mathbf{n}.\nabla u+ a(u) u=0 , & x\in\partial \Omega,\\\mathbf{n}.\nabla v+ b(v) v=0 , & x\in\partial \Omega,\\\end{array}\right.$$where $\Omega$ is a bounded smooth domain of $\mathbb{R}^{N},$ $\Delta$ is the Laplacian operator, $\lambda>0$ is a parameter, $f,g$ belong to a class of non-negative functions that have a combined sublinear effect at $\infty,$ and $a,b: [0,\infty) \rightarrow (0,\infty)$ are nondecreasing $C^{1}$ functions. We establish our existence and multiplicity results by the method of sub-- and supersolutions.