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 8kpc is ~1-3x10^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 ~4kpc radius (equivalent to ~0.5R_25_), while the MRI energy is always smaller than the turbulent energy within 8kpc 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 7kpc radius, and the sum of SNe and MRI energies is able to maintain turbulence out to 8kpc radius. Due to lack of constraint in the mass accretion rate through the disk of M33, we cannot 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.