Classical wavelet theory is based on scaling (refinement) equation whose solutions are called scaling functions. A multiwavelet of multiplicity r consists of a multiscaling function and a multiwavelet function. Multiwavelets have advantages over scalar wavelets as simultaneous inclusion of properties like short support, high approximation order, symmetry and orthogonality are possible in multiwavelets. Construction of multiwavelets deals with finding a solution to the multiscaling equation. An important observation is that a multiscaling equation in the Fourier domain can be represented in terms of a matrix polynomial H(z) called Symbol function. A matrix polynomial is a matrix whose elements are univariate or multivariate polynomials. Hence the existence and properties of a multiscaling function can be studied in terms of the spectral properties of the corresponding symbol function. The major objective of the present proposal is to explore the relationship between matrix polynomials and corresponding mutiscaling functions of multiwavelets and thereby proposing an algorithm for the construction of multiwavelets satisfying desirable properties. Development of a method to construct symbol functions corresponding to compactly supported, symmetric and orthogonal multiwavelets are very crucial in various areas of science and engineering. This motivates the development of a method to construct compactly supported, symmetric and orthogonal multiscaling functions and the corresponding multiwavelets using standard pairs. Also, we aim to construct multiscaling functions which are component wise symmetric rather than a single point. Other desirable properties of wavelets like high approximation order, high smoothness, etc. can also be studied from the perspective of matrix polynomial theory. Extra conditions on a standard pair so that the multiscaling function acquire these properties also, have to be derived. Literature includes the construction of pseudo biorthogonal pair of multiscaling functions which are constructed using cascade algorithm. However, for the convergence of the cascade algorithm the coefficients of the symbol function must satisfy sum rules of order 1. This demands the need for finding conditions on the standard pair (U, T) such that the corresponding symbol function satisfies sum rules of order 1. For the convergence of cascade algorithm, the joint spectral radius must be less than 1. Even if we construct a symbol function corresponding to a multiscaling function, the joint spectral radius need not be less than 1. We have to find the extra conditions on a standard pair (U, T) such that the joint spectral radius is less than 1. The conditions on a standard pair such that the transition matrix satisfies the condition E (1) should also be derived, since this ensures the constructed multiscaling function is square integrable.