Ordinary and partial differential equations are extensively used in different branches of science and engineering to model wide variety of phenomena, such as diffusion, stability, wave propagation, population growth and chemical reactions, to mention just a few. For most practical problems these differential equations cannot be solved analytically and numerical techniques must be employed. This paper develops methodology for applying mutli-task Gaussian processes to numerical solution of ordinary (ODEs) and partial (PDEs) differential equations. For different classes of ODEs and PDEs a detailed evaluation of the accuracy of the proposed methodology is presented by comparing the obtained numerical solutions with the corresponding exact analytical ones.