An algorithm is presented for constructing orthogonal multiscaling functions and multiwavelets with higher approximation order in terms of any given orthogonal multiscaling functions. That is, let $\Phi (x) = [\phi _{1}(x), \phi _{2}(x),\ldots , \phi _{r}(x)]^{T} \in (L^{2}(R))^{r}$ be an orthogonal multiscaling function with multiplicity $r$ and approximation order $m$. We can construct a new orthogonal multiscaling function $\Phi ^{new}(x) = [\Phi ^{T} (x), \phi _{r+1}(x), \phi _{r+2}(x),\ldots ,\phi _{r+s}(x)]^{T}$ with approximation order $n(n > m)$. Namely, we raise approximation order of a given multiscaling function by increasing its multiplicity. Corresponding to the new orthogonal multiscaling function $\Phi ^{new}(x)$, orthogonal multiwavelet $\Psi ^{new}(x)$ is constructed. In particular, the spacial case that $r = s$ is discussed. Finally, we give an example illustrating how to use our method to construct an orthogonal multiscaling function with higher approximation order and its corresponding multiwavelet.