I would like to create a 3D model inspired by Maxwell's thermodynamic surface, which shows the energy as function of entropy and volume for some "fictitious substance with water-like properties".
I would like to plot $U(S,V)$ for Van der Waals gas, which requires a change of variables from the known equation of state.
What is $U(S,V)$ for Van der Waals gas? (pre-Maxwell construction)
You require something more.
The known equation of state provides a relationship between pressure, temperature, and molar volume $(v)$: $$ \frac{p}{T}=\frac{R}{v-b} -\frac{a}{v^2} \frac{1}{T} $$ From that equation alone, it is impossible to get a unique fundamental function, either the molar entropy, $s(u,v)$ or the molar energy, $u(s,v)$. The reason is that temperature and pressure are connected to the first partial derivatives of one of these fundamental equations. The resulting equation will determine the molar entropy only within an arbitrary energy function.
Therefore, in addition to the usual van der Waals equation of state, one has to provide a thermal equation. The simplest is the relation between internal energy, volume, and temperature. In Callen's textbook (I strongly recommend it to get a better understanding of the formal structure of Thermodynamics), there is a worked example where that molar energy is assumed to be connected to the temperature and volume through the following relation $$ \frac{1}{T}=\frac{cR}{u+a/v}, $$ where $c$ is a constant.
From these two equations and using the Gibbs-Duhem equation, one gets the fundamental equation for the molar entropy: $$ s = R\log\left[ \left(v-b\right) \left( u+\frac{a}{v} \right)^c \right] + s_0. \tag{1} $$ where $s_0$ is an arbitrary additive constant.
It is trivial to invert equation $(1)$ to get $u(s,v)$.
However, as expected, this is not a convex function over the whole domain ( $v>b$, and $s>0$). There is a convex intruder that an affine region (the coexistence region) should substitute via the usual Maxwell construction or a double Legendre-Fenchel transformation.
If you only want the gas region, i.e., no vapor/liquid equilibrium, you can calculate $U(S,V)$ as follows:
If you substitute the van der Waals pressure in the equations for $S$, $H$ and $U$, you will obtain these properties as a function of $T$ and $V$. Basically, you have the result in parametric form with $T$ as the parameter.
$T_0$ and $P_0$ define the reference state, that's where $S=0$ and $H=0$. You can pick them arbitrarily, the choice will simply shift the $U(S,V)$ surface up and down by a fixed amount.
For simplicity $C_P$ is assumed constant. If you want to include its dependence in temperature, make the substitutions $$ C_P\ln\frac{T}{T_0} \to \int_{T_0}^T\frac{C_P}{T}dT,\quad C_P(T_0-T)\to \int_{T_0}^T C_P dT $$
Numerical procedure
This will give you a table of $U$ as a function of $U$ and $V$. This will also work in the liquid region, but since the van der Waals equation does not work in the two-phase region, $V_L<V<V_V$, you would need to determine $V_L$ and $V_V$ as a function of $P^\text{sat}$, which is a bit more involved (but doable).
