An algorithm for the solution of nonlinear systems of parabolic partial differential equations is presented. In this algorithm the second-order spatial derivatives are approximated by a three-point, five-point or seven-point central differencing scheme. The first-order spatial derivatives are approximated by hybrid schemes that combine central differencing with upwinding schemes that are based on two, three, four and five points interpolating polynomial formulas. The system of nonlinear ODEs that result from the spatial discretization is solved using a variable time step/variable order algorithm based on the backward differentiation formulas. Due to the inclusion of higher-order spatial differencing and upwinding, the resultant algorithm is capable of solving a wide class of nonlinear systems of one-dimensional parabolic PDEs efficiently. Realistic numerical examples are given with emphasis on the convection diffusion equation.