Partial Differential Equations (PDEs) are fundamental mathematical models for describing physical phenomena, yet most PDEs of practical interest cannot be solved analytically and require numerical approximations. The feasibility of such numerical methods, however, is ultimately constrained by the limitations of existing computational models. Since digital computers constitute the primary physical realizations of numerical computations, and Turing machines define their theoretical limits, the question of Turing computability of PDE solutions arises as a problem of fundamental theoretical significance. The Turing computability of PDE solutions provides a rigorous framework to distinguish equations that are, in principle, effectively solvable from those that inherently encode undecidable or non-computable behavior. Once computability has been established, complexity theory extends the analysis by quantifying the computational resources required to approximate the corresponding PDE solutions. In this work, we present a novel framework based on least-squares variational formulations and their associated gradient flows to analyze the computability and complexity of PDE solutions from an optimization perspective. Our approach enables the approximation of PDE solution operators via discrete gradient flows, linking structural properties of the PDE, such as coercivity, ellipticity, and convexity, to the complexity of their solutions. Within this setting, we characterize representation- and discretization-dependent sufficient conditions for regimes where PDEs admit effective, polynomial-time numerical approximations, as well as those exhibiting complexity blowup, where the input data possess polynomial-time complexity, yet the solution itself scales super-polynomially. In summary, this paper develops a general variational framework for analyzing the computability and computational complexity of classes of PDE solutions. The results provide insights into how structural properties of the underlying PDE and the regularity of the solutions influence its complexity, by establishing sufficient conditions for computability and complexity bounds. Beyond the theoretical characterization, the framework provides principled guidelines for the design of effective numerical methods and contributes to the understanding of the fundamental limitations of digital computation for PDE problems.