It seems to me that you have already done most of the work, writing the problem as an optimization problem for multiple variable. In general calling $\theta$ the actual value of the area as a function of the different points chosen.
$\theta(x_1,x_2,...,x_N) = x_1\cdot f(x_1) + \sum_{I=2}^N{f(x_I)\cdot(x_I-x_{I-1})}$
And we want to find the maximum of theta, so we search the point with all partial derivatives equal to zero
$\begin{cases}\frac{\partial \theta }{\partial x_1} = f(x_1) + x_1 \cdot f^{'}(x_1) - 0 \cdot f^{'}(x_1) - f(x_2) \\ \frac{\partial \theta }{\partial x_2} = f(x_2) + x_2 \cdot f^{'}(x_2) - x_1 \cdot f^{'}(x_2) - f(x_3) \\.... \\ \frac{\partial \theta }{\partial x_N} = f(x_N) + x_N \cdot f^{'}(x_N) - x_{N-1} \cdot f^{'}(x_N) - 0 \end{cases}$
This in addition to the conditions
$ \begin{cases}
x_J>x_I & \text{for $J>I$} \\
x_I > 0 & \forall I
\end{cases}
$
will give you a system of equation that you can solve for the value of $x_I$. The system may not always be solvable, depending on the definition of $f$. Also you should check that the chosen point is a maxima and not a minima.
Some examples:
For $f= 1-x$ the solution is always to split in equal part.
For $f= 1-x^2$ and $N=2$ the two points are
$ \begin{cases}
x_1 = \sqrt{\frac{1}{9-2\sqrt3}} \approx 0.425\\
x_2 = \sqrt{\frac{9}{9-2\sqrt3}} \approx 0.736
\end{cases}$