Many fields of science and engineering have used wavelet functions. They are established from expansion of a single mother wavelet function. Boubaker wavelet functions are presented in this paper based on the important properties of Boubaker polynomials. The research goal of this article is to drive a Boubaker wavelets operation matrix of integration in general formulas. Then an approximate solution method for solving a singular initial value problem is presented using Boubaker wavelets along the obtained operational matrix of integration. The importance of this method is that it converts a singular initial value problem in order to solve algebraic examples as a system. The process is based on reducing by means of integration the original problem into integral equations using a Boubaker wavelets operation matrix of integration to predict the integral equation. Illustrative experiments are included. In addition, computational results obtained by a Boubaker wavelets operation matrix of integration are compared with the exact solutions and other existing methods.