This paper introduces a novel approach for solving systems of boundary value problems (BVPs) by employing the recently developed Discontinuous Galerkin (DG) method, which removes the necessity for auxiliary variables. This marks the initial installment in a sequence of publications dedicated to exploring DG methods for solving partial differential equations (PDEs). In fact, through a systematic application of the DG method to each spatial variable within the PDE, employing the method of lines, we convert the initial problem into a system of ordinary differential equations (ODEs). In the current study, we developed a global error analysis of the DG method applied to systems of ODEs. Our analysis shows that using p-degree piecewise polynomials and h-mesh step size, the DG solutions achieve optimal O(h p+1 ) convergence rates in the L 2 -norm.