  We report on high-order harmonic generation in the two-band superconductor
MgB$_2$ driven by intense terahertz electromagnetic pulses. Third- and
fifth-order harmonics are resolved in time domain and investigated as a
function of temperature and in applied magnetic fields crossing the
superconducting phase boundary. The high-order harmonics in the superconducting
phase reflects nonequilibrium dynamics of the superconducting order parameter
in MgB$_2$, which is probed via nonlinear coupling to the terahertz field. The
observed temperature and field dependence of the nonlinear response allows to
establish the superconducting phase diagram.
