  In this paper we introduce a multilevel Picard approximation algorithm for
general semilinear parabolic PDEs with gradient-dependent nonlinearities whose
coefficient functions do not need to be constant. We also provide a full
convergence and complexity analysis of our algorithm. To obtain our main
results, we consider a particular stochastic fixed-point equation (SFPE)
motivated by the Feynman-Kac representation and the Bismut-Elworthy-Li formula.
We show that the PDE under consideration has a unique viscosity solution which
coincides with the first component of the unique solution of the stochastic
fixed-point equation. Moreover, the gradient of the unique viscosity solution
of the PDE exists and coincides with the second component of the unique
solution of the stochastic fixed-point equation. Furthermore, we also provide a
numerical example in up to $300$ dimensions to demonstrate the practical
applicability of our multilevel Picard algorithm.
