As far as I understand it, if you do not use the NormControl, the time step is adapted so that the maximum error on any of the solution components is below the tolerance threshold.
If, on the other hand, you use the NormControl option, then the time step is adapted so that the overall error norm is lower than the tolerance threshold. This is slightly less stringent as the biggest soltion component error is "miwed" with the other errors. Let's take a worst-case scenario: you have $N-1$ variables that follow an ODE of the type $y' = a$ with $a$ a constant. Runge-Kutta methods integrate this exactly so the error on these components is zero to machine precision. Now if you had one solution component (the $N$-th variable) which follows a nonlinear ODE, the integration error estimate $e_n$ on this variable will be nonzero.
Without the NormControl option, the time step will be adjusted roughly as follows:
$$ \Delta t \approx \left( \dfrac{|e_N|}{atol + rtol |y_N|} \right)^{\frac{1}{q+1}}$$ with $q$ the order of the integration error estimate, which is most likely equal to $p-1$, with $p$ the order of the method.
With the NormControl option set, it will be adjusted as:
$$ \Delta t \approx \left( \dfrac{ \sqrt{\sum\limits_{i=1}^{N} |e_i|^2}}{atol + rtol \sqrt{\sum\limits_{i=1}^{N} |y_i|^2}} \right)^{\frac{1}{q+1}}$$
which, in our worst-case scenario, would be:
$$ \Delta t \approx \left( \dfrac{ |e_n|}{atol + rtol \sqrt{\sum\limits_{i=1}^{N} |y_i|^2}} \right)^{\frac{1}{q+1}}$$
If we compare this to the very first equation of this answer, we see that the numerator is the same, but the denominator is larger, as $||y|| \geq |y_n|$. So the integration error will be perceived as lower with NormControl, hence the time step will be larger and the simulation quicker. But this also means that the solution quality is slightly worse.
So in your figure, I think the solution without NormControl is the better one in terms of solution exactitude. You could check that by lowering the tolerance down to very fine levels (say 1e-12) to see what the "exact" solution is. Also, you can compare how the time step evolves between the two cases to get more insight on this behaviour.
I am not quite sure, but I think the error control based on the norm of the error may be more gentle with the time step variations, whereas the solution based on the maximum error component may have bigger "jumps" in the error estimate as the solution evolves.
Also, I usually use the norm $||y|| = ||y||_{Matlab} / \sqrt{N}$, i.e. a root mean square error. This way, for discretized PDEs, the error does not "fictiously" grow as the number of mesh points increases. However, single large error components may not be represented well as they will be averaged with the lower components.
The norm in Matlab does not seem to have the factor $1/\sqrt{N}$.
Therefore, this should mean that, if your system is sufficiently large, the NormControl option will actually be more stringent on the time step. Indeed, say we duplicate a scalar ODE $N$ times, then the numerator of the above time step formula will be:
$\sqrt{\sum\limits_{i=1}^{N} |e_i|^2} = \sqrt{N |e_1|^2} = \sqrt{N} |e_1|$ which will be higher than $max(|e|)$. Thus I am guessing that your test case had a number of variables $N$ reasonably low.
In any case, the choice of using or not using this NormControl option should not affect wildly your result, otherwise this means that your integration tolerances are too large.
