  In the DARKexp framework for collisionless isotropic relaxation of
self--gravitating matter, the central object is the differential energy
distribution $n(E)$, which takes a maximum--entropy form proportional to
$\exp[-\beta(E - \Phi(0))] - 1$, $\Phi(0)$ being the depth of the potential
well and $\beta$ the standard Lagrange multiplier. Then the first and quite
non--trivial problem consists in the determination of an ergodic phase--space
distribution which reproduces this $n(E)$. In this work we present a very
extensive and accurate numerical solution of such DARKexp problem for systems
with cored mass density and finite size. This solution holds throughout the
energy interval $\Phi(0)\le E\le 0$ and is double--valued for a certain
interval of $\beta$. The size of the system represents a unique identifier for
each member of this solution family and diverges as $\beta$ approaches a
specific value. In this limit, the tail of the mass density $\rho(r)$ dies off
as $r^{-4}$, while at small radii it always starts off linearly in $r$, that is
$\rho(r)-\rho(0)\propto r$.
