In this paper, we consider theoretical and numerical properties of a nonlinear boundary-value problem which is strongly related to the well-known Gelfand–Bratu model with parameter λ. When approximating the nonlinear term in the model via a Taylor expansion, we are able to find new types of solutions and multiplicities, depending on the final index N in the expansion. The number of solutions may vary from 0, 1, 2 to ∞. In the latter case of infinitely many solutions, we find both periodic and semi-periodic solutions. Numerical experiments using a non-standard finite-difference (NSFD) approximation illustrate all these aspects. We also show the difference in accuracy for different denominator functions in NSFD when applied to this model. A full classification is given of all possible cases depending on the parameters N and λ.