A scaling function is the solution to a dilation equation Φ(t) = ΣckΦ(2t - k), in which the coefficients come from a low-pass filter. The coefficients in the wavelet W(t) = ΣdkΦ(2t - k) come from a high-pass filter. When these coefficients are matrices, Φ and W are vectors: there are two or more scaling functions and an equal number of wavelets. By dilation and translation of the wavelets, we have an orthogonal basis Wijk = Wi(2jt - k) for all functions of finite energy. These "multiwavelets" open new possibilities. They can be shorter, with more vanishing moments, than single wavelets. They can be symmetric, which is impossible for scalar wavelets (except for Haar's). We determine the conditions to impose on the matrix coefficients ck in the design of multiwavelets, and we construct a new pair of piecewise linear orthogonal wavelets with two vanishing moments.