A key and expensive part of coupled atmospheric chemistry-climate model simulations is the integration of gas phase chemistry, which involves dozens of species and hundreds of reactions. These species and reactions form a highly-coupled network of Differential Equations (DEs). There exists orders of magnitude variability in the lifetimes of the different species present in the atmosphere and so solving these DEs to obtain robust numerical solutions poses a “stiff problem”. With newer models having more species and increased complexity it is now becoming increasingly important to have chemistry solving schemes that reduce time but maintain accuracy. While a sound way to handle stiff systems is by using implicit DE solvers, the computational costs for such solvers are high due to internal iterative algorithms (e.g., Newton-Raphson methods). Here we propose an approach for implicit DE solvers that improves their convergence speed and robustness with relatively small modification in the code. We achieve this by blending the existing Newton-Raphson (NR) method with Quasi-Newton (QN) methods, whereby the QN routine is called only on selected iterations of the solver. We test our approach with numerical experiments on the UK Chemistry and Aerosol (UKCA) model, part of the UK Met Office Unified Model suite, run in both an idealized box-model environment and under realistic 3D atmospheric conditions. The box model tests reveal that the proposed method reduces the time spent in the solver routines significantly, with each QN call costing 27% of a call to the full NR routine. A series of experiments over a range of chemical environments was conducted with the box-model to find the optimal iteration steps to call the QN routine which result in the greatest reduction in the total number of NR iterations whilst minimising the chance of causing instabilities and maintaining solver accuracy. The 3D simulations show that our moderate modification, by means of using a blended method for the chemistry solver, speeds up the chemistry routines by around 13%, resulting in a net improvement in overall run-time of the full model by approximately 3 % with negligible loss in the accuracy. The blended QN method also improves the robustness of the solver, reducing the number of grid cells which fail to converge after 50 iterations by 40%. The relative differences in chemical concentrations between the control run and that using the blended QN method are of order ~10−7 for longer lived species, such as ozone, and below the threshold for solver convergence (10-4) almost everywhere for shorter lived species such as the hydroxyl radical.