A hybrid efficient and highly accurate spectral matrix technique is adapted for numerical treatments of a class of two-pint boundary value problems (BVPs) with singularity and strong nonlinearity. The underlying model is a reaction-diffusion equation arising in the modeling of biomedical, chemical, and physical applications based on the assumptions of Michaelis–Menten kinetics for enzymatic reactions. The manuscript presents a highly computational spectral collocation algorithm for the model combined with the quasilinearization method (QLM) to make the proposed technique more efficient than the corresponding direct spectral collocation algorithm. A novel class of polynomials introduced by S.K. Chatterjea is used in the spectral method. A detailed proof of the convergence analysis of the Chatterjea polynomials (ChPs) is given in the L2 norm. Different numerical examples for substrate concentrations with all values of parameters are performed for the case of planar and spherical shapes of enzymes. For validation, these results are compared with those obtained via wavelet-based procedures and the Adomian decomposition scheme. To further improve the approximate solutions obtained by the QLM–ChPs method, the technique of error correction is introduced and applied based on the concept of residual error function. Overall, the presented results with exponential convergence indicate that the QLM–ChPs algorithm is simple and flexible enough to be applicable in solving many similar problems in science and engineering.