In this thesis, we study the limiting distributions of linear systems of partial differential equations with subordinated Gaussian random initial data. When the initial data is non-random, the solutions of the linear systems are given by the convolution of the Green kernels and the initial data. Therefore, the evolution of the solutions is totally determined by their initial data. However, the information regarding the initial data is obtained through some process of measurement, resulting in measurement error. In our work, we use the second-order homogeneous random field to model these measurement error and apply the spectral representation method to study the covariance matrix functions of the random solution vector fields. In view of that the solution fields can be thought of as the weighted sum of correlated random variables, we will also consider the limiting distributions of the random solution fields from di↵erent viewpoints, including macroscopic scales and microscopic scales. When the random initial data is weakly dependent, our results can be thought of as a generalized central limit theorem. There are two contributions for the new results. The first one is that the initial data is modeled by two cross-correlated subordinated Gaussian random fields. We use the method of Feynman diagrams to analyze the asymptotic behavior of the covariance matrix function of the random solution field induced by the random initial data. Second, the limit of the random solution vector field under the macroscopic/microscopic coordinate systems is represented by a L2-convergent series of mutually independent Gaussian random fields. We also study the limiting distributions of the solution vector field when its random initial data is long-range dependent. Compared to the previous case, the limiting law of the rescaled solution vector field is non-Gaussian, which is represented by multiple Wiener integrals. In contrast to the existing mathematical literature we found that there is a competition relationship between the effect coming from two components of the random initial data. That is, one of the two components of the random initial data will be determined dominantly the structure of the limiting distribution of the random part of the solution vector field.