In this article, we address the design and analysis of multivariable Newton-based extremum seeking (ES) for static locally quadratic maps subject to actuation dynamics governed by diffusion partial differential equations (PDEs). Multi-input systems with distinct diffusion coefficients in each individual input channel are dealt with. The phase compensation of the dither signals is handled as a trajectory generation problem and the inclusion of a multivariable diffusion feedback controller with a perturbationbased (averaging-based) estimate of the Hessian's inverse allow to obtain local exponential convergence results to a small neighborhood of the optimal point. The stability analysis is carried out using backstepping transformation and averaging in infinite dimensions, capturing the infinite-dimensional state due the diffusion PDEs. In addition, the generalization of the results for different classes of parabolic (reaction-advection-diffusion) PDEs, wave equations and/or first-order hyperbolic (transport-dominated) PDEs is also discussed. As for ordinary differential equations case, the proposed Newton approach removes the dependence of the algorithm's convergence rate on the unknown Hessian of the nonlinear map to be optimized, being user-assignable unlike the gradient algorithm. A numerical example illustrates the performance of our multivariable ES for compensating PDE-based systems.