In this thesis, we will investigate McKean-Vlasov SDEs (McKV-SDEs) on Rd,: ; where coefficient functions b and σ satisfy sufficient regularity conditions and fWtgt2[0;T ] is a Wiener process. These SDEs correspond to a class of deterministic non-local PDEs. The principal aim of the first part is to present Multilevel Monte Carlo (MLMC) schemes for the McKV-SDEs. To overcome challenges due to the dependence of coefficients on the measure, we work with Picard iteration. There are two different ways to proceed. The first way is to address the McKV-SDEs with interacting kernels directly by combining MLMC and Picard. The MLMC is used to represent the empirical densities of the mean-fields at each Picard step. This iterative MLMC approach reduces the computational complexity of calculating expectations by an order of magnitude. However, we can also link the McKV-SDEs with interacting kernels to that with non-interacting kernels by projection and then iteratively solve the simpler by MLMC method. In each Picard iteration, the MLMC estimator can approximate a few mean-fields directly. This iterative MLMC approach via projection reduces the computational complexity of calculating expectations hugely by three orders of magnitude. In the second part, the main purpose is to demonstrate the plausibility of applying deep learning technique to several types of random linear PDE systems. We design learning algorithms by using probabilistic representation and Feynman- Kac formula is crucial for deriving the recursive relationships on constructing the loss function in each training session.