As long as long-range interactions (e.g., gravity) are negligible, the internal energy $U$ of a homogeneous system is an extensive quantity. In such a case, the second law of thermodynamics in fact guarantees that $U$ is a convex function of other extensive quantities it depends on, e.g., $S$, $V$, and $N$.
An important observation is that given a thermodynamic system, a tiny part of it should be in equilibrium with the rest. Suppose that we take this tiny part to be negligibly small compared to the entire system, and simply call it "the system." The rest of the entire system acts as a reservoir with some fixed temperature $T_{0}$, pressure $P_{0}$, and chemical potential $\mu_{0}$, so let's call this "the reservoir." When the system and the reservoir are in thermodynamic equilibrium, the second law of thermodynamics demands that the total entropy be maximized. Notice that the change in the total entropy is given by
\begin{equation}
\Delta S_{\mathrm{total}} = \Delta S + \frac{Q_{\mathrm{reservoir}}}{T_{0}} = \Delta S - \frac{\Delta U + P_{0}\Delta V - \mu_{0}\Delta N}{T_{0}} = -\frac{1}{T_{0}} \Delta \Omega,
\end{equation}
where $\Omega = U-T_{0}S + P_{0}V-\mu_{0}N$. Hence, maximizing $S_{\mathrm{total}}$ is equivalent to minimizing $\Omega$. Here, $\Omega$ has three, as opposed to four, independent variables because there is a relation between $U$, $S$, $V$, and $N$. Taking the independent parameters to be $S$, $V$, and $N$, we have
\begin{equation}
\Omega(S,V,N) = U(S,V,N) - T_{0}S + P_{0} V - \mu_{0} N.
\end{equation}
A consistency condition is that for an arbitrary $T_{0}$, $P_{0}$, and $\mu_{0}$, there exists a unique minimum of $\Omega(S,V,N)$ with respect to its independent variables, viz., the system should be able to equilibrate with the reservoir. It then follows that the function $U(S,V,N)$ is convex.
So far, $S$, $V$, $N$, and $U(S,V,N)$ are properties of a small system that is a part of a much larger body. However, for a homogeneous system, the internal energy of the entire body must have the same form as any small part of it. Therefore, what I have shown equally holds for an arbitrary homogeneous system.
[Taking the independent variables of $\Omega$ to be $U$, $V$, and $N$ would lead to the concavity of the function $S(U,V,N)$.]
