In this paper, we will introduce an algorithm for obtaining integrals of the form \begin{equation*} \int_0^x t^m \varphi(t)dt,\quad m\in \mathbb{N}\cup\{0\}, \end{equation*} where $\varphi$ is the scaling functions of Daubechies wavelet. In order to obtain these integrals in dyadic points for $x$'s, we have to solve a linear system. We show that these matrices which obtaining by integrals involving Daubechies scaling functions and plynomials are bounded. Also, we will investigate, sparseness, well-conditioning and strictly diagonal dominant of matrices of these systems.
