  We compute the quasinormal modes of massive scalar and Dirac fields within
the framework of asymptotically de Sitter black holes in Euler-Heisenberg
non-linear electrodynamics. We pay particular attention to the regime $\mu
M/m_{P}^2 \gg 1$, where $\mu$ and $M$ denote the masses of the field and the
black hole, respectively, and $m_{P}$ represents the Planck mass, covering a
range from primordial to large astrophysical black holes. Through time-domain
integration, we demonstrate that, contrary to the asymptotically flat case, the
quasinormal modes also dictate the asymptotic decay of fields. Employing the
6th order WKB formula, we derive a precise analytic approximation for
quasinormal modes in the regime $\mu M/m_{P}^2 \gg 1$ without resorting to
expansion in terms of the inverse multipole number. This analytic expression
takes on a concise form in the limit of linear electrodynamics, represented by
the Reissner-Nordstr\"om black holes.
