We utilize the multi-wavelength data of M33 to study the origin of turbulence in its interstellar medium. We find that the HI turbulent energy surface density inside 8 kpc is ∼ 1 − 3 × 10 46 erg pc −2 , and has no strong dependence on galactocentric radius because of the lack of variation in HI surface density and HI velocity dispersion. Then, we consider the energies injected by supernovae (SNe), the magneto-rotational instability (MRI), and the gravity-driven turbulence from accreted materials as the sources of turbulent energy. For a constant dissipation time of turbulence, the SNe energy can maintain turbulence inside ∼ 4 kpc radius (equivalent to ∼ 0.5 R 25 ), while the MRI energy is always smaller than the turbulent energy within 8 kpc radius. However, when we let the dissipation time to be equal to the crossing time of turbulence across the HI scale-height, the SNe energy is enough to maintain turbulence out to 7 kpc radius, and the sum of SNe and MRI energies is able to maintain turbulence out to 8 kpc radius. Due to lack of constraint in the mass accretion rate through the disk of M33, we can not rule out the accretion driven turbulence as a possible source of energy. Furthermore, by resolving individual Giant Molecular Clouds in M33, we also show that the SNe energy can maintain turbulence within individual molecular clouds with ∼ 1% of coupling efficiency. This result strengthens the proposition that stellar feedback is an important source of energy to maintain turbulence in nearby galaxies.