  The aim of this paper is to study the time stepping scheme for approximately
solving the subdiffusion equation with a weakly singular source term. In this
case, many popular time stepping schemes, including the correction of
high-order BDF methods, may lose their high-order accuracy. To fill in this
gap, in this paper, we develop a novel time stepping scheme, where the source
term is regularized by using a $k$-fold integral-derivative and the equation is
discretized by using a modified BDF2 convolution quadrature. We prove that the
proposed time stepping scheme is second-order, even if the source term is
nonsmooth in time and incompatible with the initial data. Numerical results are
presented to support the theoretical results.
