A formula for integrals (primitives) of Bernstein polynomials in terms of sums of Bernstein polynomials is given. For a very special class of
th order linear ordinary differential equations, subject to Dirichlet boundary conditions, a numerical approximation method is proposed by integrating the differential equation
-times and using Bernstein polynomials as basis functions for a Galerkin method. The above-mentioned formula for integrals of Bernstein polynomials is used to calculate the corresponding matrix elements. Some numerical examples illustrate the approach.