rather "shorter" answer would be going directly from conservation of 4-momentum and it's squares.
Let us denote 4-momenta of particles by N, P, k, k' for neutron, proton, electron and electron anti-neutrino in this respective order.
Total 4-momentum of an isolated system is conserved, thusly it holds
$N - P - k - k' = 0$, or $(N - k) = (P + k')$.
Note that when either LHS or RHS are squared (i.e. contracted with metric) you get and invariant. A quantity which is conserved under Lorentzian transform is precisely the same in value in every inertial frame.
1] Let's square both LHS and RHS.
2] I am using $g_{ab} = diag(1, -1, -1, -1)$ and thus 4-momentum is normalized to $+m^2$
3] We wish to get, after squaring, on RHS just linear term in $E_e$ (electron energy). It is convenient to use neutron rest frame, where in particular $N = \left(E_n, \overrightarrow{0}\right) = \left(m_n, \overrightarrow{0}\right)$
$$RHS^2 = N^2 -2N.k + k^2 = m_n^2 + m_e^2 - 2m_n E_e$$
$$LHS^2 = P^2 + k^2 + 2P.k = m_p^2 + m_\nu^2 + 2P.k \ge m_p^2 + m_\nu^2 $$
Together ($LHS^2 = RHS^2$ due to invariance) that is
$$ -2m_n E_e \ge - m_n^2 - m_e^2 + m_p^2 + m_\nu^2 $$
$$ E_e \le \frac{m_n^2 - m_p^2 + m_e^2 - m_\nu^2}{2m_n} $$
Now to arrive at your result:
1] $m_e$ is orders of magnitude smaller than $m_n$ in the denominator, $m_\nu = 0$ for all practical purposes of this level of calculation. We neglect both squares of those masses.
2] $(m_n^2-m_p^2) = (m_n + m_p)(m_n - m_p) =: 2M(m_n - m_p)$, where I have introduced approximate mass of nucleons $M$ for sum of two nucleons (not for their difference). If you goof around which is fun you find this brings out around $1/1000$ error if I remember correctly.
$$ E_e \le \dots \approx \frac{2M(m_n - m_p)}{2m_n} \approx \frac{2M(m_n - m_p)}{2M} \approx \frac{2M(m_n - m_p)}{2M} = m_n - m_p$$
i.e. upper energy boundary for electron is about the difference of the proton-neutron mass.
Sidenote:
If I wanted to extract information about instead $E_p$, I would have started from neutron rest frame and squared $(N-P)$ to get linear term in $E_p$. I would need to use also for LHS electron-neutrino COM frame to get $\overrightarrow{k'} = -\overrightarrow{k}$ and use it in inequality of the squares similarly as above.