The method [whuber](https://stats.stackexchange.com/users/919/whuber) has used in his excellent answer is a common optimisation "trick" that involves extending the likelihood function to allow real values of $N$, and then using the concavity of the log-likelihood to show that the discrete maximising value is one of the discrete values on either side of a continuous optima.  This is one commonly used method in discrete MLE problems involving a concave log-likelihood function.  Its value lies in the fact that it is usually possible to get a simple closed-form expression for the continuous optima.

For completeness, in this answer I will show you an alternative method, which uses discrete calculus using the [forward-difference operator](https://en.wikipedia.org/wiki/Finite_difference).  The log-likelihood function for this problem is the discrete function:

$$\ell_y(N) = -\frac{1}{2} \Bigg[ \ln (2 \pi) + \ln (\sigma^2) + \ln (N) + \frac{(y-N\mu)^2}{N\sigma^2} \Bigg]
\quad \quad \quad
\text{for } N \in \mathbb{N}.$$

The first forward-difference of the log-likelihood is:

$$\begin{equation} \begin{aligned}
\Delta \ell_y(N) 
&= -\frac{1}{2} \Bigg[ \ln (N+1) - \ln (N) + \frac{(y-N\mu - \mu)^2}{(N+1)\sigma^2} - \frac{(y-N\mu)^2}{N\sigma^2} \Bigg] \\[6pt]
&= -\frac{1}{2} \Bigg[ \ln \Big( \frac{N+1}{N} \Big) + \frac{N(y-N\mu - \mu)^2 - (N+1)(y-N\mu)^2}{N(N+1)\sigma^2} \Bigg] \\[6pt]
&= -\frac{1}{2} \Bigg[ \ln \Big( \frac{N+1}{N} \Big) + \frac{[N(y-N\mu)^2 -2N(y-N\mu) \mu + N \mu^2] - [N(y-N\mu)^2 + (y-N\mu)^2]}{N(N+1)\sigma^2} \Bigg] \\[6pt]
&= -\frac{1}{2} \Bigg[ \ln \Big( \frac{N+1}{N} \Big) - \frac{(y + N \mu)(y-N\mu) - N \mu^2}{N(N+1)\sigma^2} \Bigg]. \\[6pt]
\end{aligned} \end{equation}$$

With a bit of algebra, the second forward-difference can be shown to be:

$$\begin{equation} \begin{aligned}
\Delta^2 \ell_y(N) 
&= -\frac{1}{2} \Bigg[ \ln \Big( \frac{N+2}{N} \Big) + \frac{2 N (N+1) \mu^2 + 2(y + N \mu)(y-N\mu)}{N(N+1)(N+2)\sigma^2} \Bigg] < 0. \\[6pt]
\end{aligned} \end{equation}$$

This shows that the log-likelihood function is concave, so its smallest maximising point $\hat{N}$ will be:

$$\begin{equation} \begin{aligned}
\hat{N} 
&= \min \{ N \in \mathbb{N} | \Delta \ell_y(N) \leqslant 0 \} \\[6pt]
&= \min \Big\{ N \in \mathbb{N} \Big| \ln \Big( \frac{N+1}{N} \Big) \geqslant \frac{(y + N \mu)(y-N\mu) - N \mu^2}{N(N+1)\sigma^2} \Big\}.
\end{aligned} \end{equation}$$

(The next value will also be a maximising point if and only if $\Delta \ell_y(N) = 0$.)  The MLE (either the smallest, or the whole set) can be programmed as a function via a simple ```while``` loop, and this should be able to give you the solution pretty quickly.  I will leave the programming part as an exercise.