In this paper, we 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 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 investigate, sparseness, well-conditioning and strictly diagonal dominant of matrices of these systems.