Multiwavelets are wavelets with multiplicity r, that is r scaling functions and r wavelets, which define multiresolution analysis similar to scalar wavelets. They are advantageous over scalar wavelets since they simultaneously posse symmetry and orthogonality. In this work, a new method for constructing multiwavelets with any approximation order is presented. The method involves the derivation of a matrix equation for the desired approximation order. The condition for approximation order is similar to the conditions in the scalar case. Generalized left eigenvectors give the combinations of scaling functions required to reconstruct the desired spline or super function. The method is demonstrated by constructing a specific class of symmetric and non-symmetric multiwavelets with different approximation orders, which include Geranimo-Hardin-Massopust (GHM), Daubechies and Alperts like multi-wavelets, as parameterized solutions. All multi-wavelets constructed in this work, posses the good properties of orthogonality, approximation order and short support.