This paper deals with Physics-Informed Gaussian Processes (GPs), a probabilistic method to solve linear parametric Partial Differential Equations (PDEs). In fact, we build upon \cite{RAISSI2017683} and use GP priors to infer statistical estimators to the PDE solution given possibly noisy observations. Moreover, we are particularly interested in the big data regime, in which we deal with a large number of observations $n$. Nonetheless, it is well-known that GPs suffer from computational burdens as it implies inverting a $n^2$-matrix. We then show how to leverage sketching approximation, such as the Nystr\"om approximation \cite{rudi2015less}, in order to scale Physics-Informed GPs up to big data, while guaranteeing good statistical accuracy. To the best of our knowledge, it is the first time that sketching is applied to Physics-Informed GPs or kernel methods. Finally, this paper is devoted to bridge a gap in the literature and provide convergence rates of such estimators to the PDE solution in terms of Sobolev norms \cite{teckentrup2020convergence}. Unlike \cite{Cockayne2017} that considers using the PDE's Green function as covariance kernel for the GP, or \cite{BATLLE2025113488} that assumes the stability of the differential operator and the solution map of the PDE between corresponding Sobolev spaces, we highlight sufficient assumptions on the target function, differential operator and covariance kernels to obtain such stability and convergence rates. In particular, we show that the Mat\'ern kernels fall into these conditions. We illustrate these points with different numerical experiments.