  We couple the L1 discretization for Caputo derivative in time with spectral
Galerkin method in space to devise a scheme that solves quasilinear
subdiffusion equations. Both the diffusivity and the source are allowed to be
nonlinear functions of the solution. We prove method's stability and
convergence with spectral accuracy in space. The temporal order depends on
solution's regularity in time. Further, we support our results with numerical
simulations that utilize parallelism for spatial discretization. Moreover, as a
side result we find asymptotic exact values of error constants along with their
remainders for discretizations of Caputo derivative and fractional integrals.
These constants are the smallest possible which improves the previously
established results from the literature.
