We couple the L1 discretization for Caputo derivative in time with spectral Galerkin method in space to devise a scheme that solves strongly nonlinear 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 order $2-\alpha$ in time and spectral accuracy in space. Further, we illustrate 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.