This paper presents a numerical method based on Fluctuationlessness Theorem for the solution of Ordinary Differential Equations over appropriately defined Hilbert Spaces. We focus on the linear differential equations in this work. The approximated solution is written in the form of an nth degree polynomial of the independent variable. The unknown coefficients are obtained by setting up a system of linear equations which satisfy the initial or boundary conditions and the differential equation at the grid points, which are constructed as the independent variable's matrix representation restricted to an n dimensional subspace of the Hilbert Space. An error comparison of the numerical solution and the MacLaurin series with the analytical solution is performed. The results show that the numerical solution obtained here converges to the analytical solution without using too many mesh points.