This article scrutinizes the 2-dimensional and boundary layer flow of magnetohydrodynamic Williamson fluid flowing on a stretchable surface with variable viscosity. The thermal and solutal rates are examined through the Cattaneo-Christov model with Joule heating, heat source/sink, and chemical reaction. The authors are motivated to conduct this study because of its practical and scientific significance in various processes, including polymer processing, textile industries, food industries, solar energy, biomedical science, wind turbine blades, oil spill clean-up, metal rolling, and forging. With the mentioned assumptions, the partial differential equations are achieved by using the basic governing laws, including momentum law, energy law, and concentration law. This non-linear system of equations is transmuted into ordinary differential equations by taking similarity transformations. The main novelty behind the conduction of this work is the numerical technique, namely the ‘Adams-Milne (Predictor-Corrector)’ method along with the Runge-Kutta technique on Matlab software, which has not previously been studied by any researcher in the literature. The analytical solution of the determined equations is not possible due to their highly non-linear nature; therefore the multistep numerical method namely the ‘Adams-Milne (Predictor-Corrector)’ method, along with the Runge-Kutta technique is used to determine the numerical results. The outcomes are noted due to numerous parameters for velocity, temperature, and concentration profiles. The explanation of graphical and numerical results is discussed here. The graphical impression of the Williamson parameter reveals that the velocity and temperature curves diminish for higher inputs of this parameter. The movement of fluid shows the declining behavior for the Hartmann number and viscosity parameter. The solutal and thermal findings due to Cattaneo-Christov heat and mass relaxation coefficients mark the reducing behaviour in respective field. The rise in reaction coefficient decreases the mass distribution. The analyses of comparison of results are also presented here.