We present a new method for solving regular boundary value problems for linear ordinary differential equations with constant coefficients (the case of variable coefficients can be adopted readily but is not treated here). Our approach works directly on the level of operators and does not transform the problem to a functional setting for determining the Green’s function. We proceed by representing operators as noncommutative polynomials, using as indeterminates basic operators like differentiation, integration, and boundary evaluation. The crucial step for solving the boundary value problem is to understand the desired Green’s operator as an oblique Moore–Penrose inverse. The resulting equations are then solved for that operator by using a suitable noncommutative Gröbner basis that reflects the essential interactions between basic operators. We have implemented our method as a Mathematica™ package, embedded in the TH ∃ OREM ∀ system developed in the group of Prof. Bruno Buchberger. We show some computations performed by this package.