I'll explain to you why a Young experiment with $N$ slits gives approximately equally spaced maxima with intensity $I_0 N^2$ (if $N\ge 3$ there are secodary maxima too, but they are small and they became negligible id $N$ is big). For the moment I will ignore problems linked to diffraction and I will focus only about far interference with $N$ coherent point sources, I will spend some words about diffraction in the end. Let's consider a row of $N$ coherent point sources separated by a distance $d$, as in the following figure
We can assume that rays light arriving in a point $P$ on the far screen are parallel, and that the field is the real part of (here, as often happens if we have to manage many waves, using complex notation is useful)
\begin{align}
\begin{split}
E_P &= E_0 e^{i (kr_1-\omega t)}+ E_0 e^{i (kr_2-\omega t)}+ \dots E_0 e^{i (kr_N-\omega t)} \\
& = E_0 e^{-i\omega t}e^{i kr_1} [1+ e^{ik(r_2-r_1)} +e^{ik(r_3-r_1)} +\dots +e^{ik(r_N-r_1)} ]
\end{split}
\end{align}
where $r_i$ is path from $i^{\textrm{th}}$ source to $P$.
If $r_1$ and $r_2$ are the two upper rays, the path difference is $r_2-r_1=d\sin \theta$ (being $\theta$ the angular position of the point $P$ on the far screen), we have also $r_3-r_1=2d\sin \theta$ and so on. By defining $\delta = kd \sin \theta$
we have then $k (r_2 - r_1) = \delta$, $k (r_3 - r_1) = 2 \delta$ and so on. So we can write
\begin{equation}\label{9331}
E_P = E_0 e^{-i\omega t} e^{ikr_1} [(e^{i\delta})^0 + (e^{i\delta})^1 + (e^{i\delta})^2 + \dots (e^{i\delta})^{N-1}]
\end{equation}
But $\sum_{n=0}^{N-1} x^n = \frac{x^N-1}{x-1}$
so what it is in square brackets is
\begin{equation}\label{7152}
\frac{e^{i\delta N}-1}{e^{i\delta}-1}
\end{equation}
It is easy to see that $\frac{ e^{\frac{iN\delta}{2}} \left[ e^{\frac{iN\delta}{2}} - e^{-\frac{iN\delta}{2}} \right] }{ e^{\frac{i\delta}{2}} \left[ e^{\frac{i\delta}{2}} - e^{-\frac{i\delta}{2}} \right]}$ is equivalent to the above expression. But this (apparently ugly) way of writing is useful because we can exploit $\sin x = \frac{e^{ix}+e^{-ix}}{2i}$ (by the way, this amazing formula, written by Euler in the middle of XVIII century, can be easily derived if we assume true Euler's formula $e^{ix}=\cos x + i \sin x$: start from generic complex number $z=a+ib$ and check that its imaginary part is $\textrm{Im}(z) = \frac{z-z^*}{2i}$, on the other hand if $z=e^{ix}$ we have $\textrm{Im}(z)=\sin x$, from this the formula). We write the expression in this way
\begin{equation}
e^{\frac{i(N-1)\delta}{2}} \frac{\sin \frac{N \delta}{2} }{ \sin \frac{\delta}{2} }
\end{equation}
Now we wrote the square bracket in a more interesting way, we can go back and say that to total electric field in $P$ is the real part of
\begin{equation}\label{9330}
E_P = E_0 e^{-i\omega t} e^{i\left[ k r_1 + \frac{(N-1)\delta}{2} \right]}
\frac{\sin \frac{N \delta}{2} }{ \sin \frac{\delta}{2} }
\end{equation}
Now, if $r_1$ is the distance of $P$ (point of the screen) from the first point source, the distance $R$ of $P$ from the center of the row of point sources (this distance can be taken as the best representative distance of $P$ from the row of sources) is $R=r_1 + \frac{1}{2} (N-1) d \sin \theta$ (because $\frac{(N-1)d}{2}$ is the distance of the first point source from the center of the row). Taking $r_1$, substituting in the above $E_P$ expression, and substituting $\delta$ definition too, we can write the electric field in this way
\begin{equation}\label{9329}
E_P = E_0 e^{i(kR-\omega t)}
\frac{\sin \frac{N \delta}{2} }{ \sin \frac{\delta}{2} }
\end{equation}
Where $R$ is the distance of the point on the screen from the center of the row of point sources. The electric field in $P$ is the real part of this expression, and it is an harmonic wave in the form $E_{max} \cos(kR-\omega t)$ where $E_{max} = E_0 \frac{\sin \frac{N \delta}{2} }{ \sin \frac{\delta}{2} } $ is a constant (please note that the term with sines is a constant for a given point $P$ on the screen, with its angular position $\theta$; consider that wave number $k$ too is ``given'', and that we assume that $E_0$ is independent from $\theta$). The intensity is directly proportional to maximum electric field: $I =\frac{\epsilon_0 c}{2} E_{max}^2 = \frac{\epsilon_0 c}{2} E_0^2 \frac{\sin^2 \frac{N \delta}{2} }{ \sin^2 \frac{\delta}{2} }$. If we set $I_0 = \frac{\epsilon_0 c}{2} E_0^2 $, we get
\begin{equation}\label{7151}
I(\theta)=I_0 \left( \frac{\sin \left( \frac{N \pi d\sin \theta}{\lambda} \right)}{\sin \left( \frac{\pi d\sin \theta}{\lambda} \right)} \right)^2
\end{equation}
Please note that we didn't exploit the hypothesis $N \gg 1$, and exploiting $\sin (2x) = 2 \sin x \cos x$ it is easy to see that this result can be exploited to find Young formula (two sources)
\begin{equation}
I (\theta) = 4 I_0 \cos^2 \left( \frac{\pi d \sin \theta}{\lambda_0 L} \right)
\end{equation}
(obviously if we are interested only to the two slit problem, this proof is not the simpler, but it is important to show that the two results are consistent).
Now we are ready t find positions and intensities of maxima. If $\theta=0$ all waves are in phase and we must have a maximum. The intensity is
\begin{equation}\label{7145}
I_0 \lim_{x \to 0} \left( \frac{\sin (Nx)}{\sin x} \right)^2 =
I_0 \left( \lim_{x \to 0} \frac{\sin (Nx)}{\sin x} \right)^2 =
I_0 \left( \lim_{x \to 0} \frac{N \cos (Nx)}{\cos x} \right)^2 = I_0 N^2
\end{equation}
where we set $\frac{\pi d \sin \theta}{\lambda} = x$ and we exploited de l'Hopital rule to solve the limit. Are there limiting values other than $\theta \to 0$ such that $x\to 0$? Surely if $\sin \theta \to 0$ we have $x \to 0$, but the only $\theta$ with physical sense are between $-\frac{\pi}{2}$ and $\frac{\pi}{2}$, so with $\theta=0$ viewing the problem under this angle has ended the possibilities. However note what follow
We have the same indeterminate limit, which is solved in the same way with de l'Hopital, not only when $x \to 0$ but also when $x\to \pi m$ with $m$ integer (remember that $N$ too is an integer, and eventual minus sign are removed when squaring)
There are no doubt that when $\theta=0$ we have a maximum of intensity on the screen (all waves are in phase) so any other point with the same intensity $I_0 N^2$ is a maximum too
We conclude that the $\theta$ values satisfying $\frac{\pi d \sin \theta}{\lambda} = m \pi$, i.e.
\begin{equation}
d \sin \theta = m \lambda \qquad m \in \mathbb{Z}
\end{equation}
are the angular position of maxima, with intensity $I_0 N^2$. If $\theta$ is small, $\sin \theta \approx \theta$ and we see that on the screen maxima are equally spaced at distance $D \Delta \theta$ being $\Delta \theta = \frac{\lambda}{d}$ ($d$ distance between point sources) and $D$ distance of the screen.
Finally, some words about diffraction problem (this explain the plots of $\frac{I}{I_0}$ in the question). We implicitly assumed $a \ll \lambda $ ($a$ amplitude of the slit: we assumed "point" sources), so that each source illuminates the screen approximately evenly (the central diffraction maximum is large).
If this is not the case, we must assume that each "point" source is not a point, and interference effect from different point of the same source force us to assume that intensity of each source has maxima and minima depending on $\theta$. Because the row of sources is small and far from screen, we can assume that this effect is the same for all sources. The intensity profile is so "modulated" by this effect, and we can account of this simply multiplying $I(\theta)$ by the squared single slit diffraction function
$\left( \frac{\sin \alpha}{\alpha} \right)^2$ with $\alpha = \frac{\pi a\sin \theta}{\lambda}$ ($a$ amplitude of the slit): the spacing between principal maxima is the same, but in this case the only one with intensity $I_0 N^2$ is the one in the center (because $\frac{\sin x}{x} < 1$ and only if $x\to 0$ we have $ \frac{\sin x}{x} \to 1$).
