I usually find it easier to use model multipoles that are surface charges on a sphere, rather than point charges on some polyhedron's vertices. These charge densities are given in general by $$\sigma_{lm}(\theta,\phi)=N\cos(m\phi)P_l^m(\cos(\theta)).$$ Thus a monopole is constant, a dipole has $\sigma=\cos(\theta)=z/r$, a quadrupole has the form $\sigma=\cos(\phi)\cos(\theta)\sin(\theta)=xz/r^2$ or $\sigma=3\cos^2(\theta)-1=\frac{2z^2-x^2-y^2}{r^2}$, and so on. This has two important advantages:
- the field produced by these surface charge densities is has pure multipolarity, unlike the point-charge multipole combinations, which always have subleading corrections. Even better, this happens outside the sphere (giving a far-field multipole) and outside it (giving a near-field multipole).
- The "and so on" actually does hold. For any $l, m$ combination you just plug the indices in and see what the charge density looks like. Even better, this covers all possible multipole configurations (whereas the point-line-square-cube game misses an important quadrupole configuration, given by dmckee in the comments, and a number of different octupole ones.)
Because of all that, let me address your question: what do hexadecapoles look like? Here $l=4$, so you have five different configurations, for $m$ from 0 through 4. In that order, they look like
$m=0$,
(which you can emulate with five charges of the appropriate magnitudes on a line),
$m=1$,
(which you can emulate with two lines of four charges),
$m=2$,
(which you can emulate with a square prism and a square),
$m=3$,
(which you can emulate with two hexagons with alternating charges), and
$m=4$,
(which you can emulate with an octagon).
For sketching purposes, it is usually enough to know that $l$ tells you how many nodal lines the charge density can have (check above!). The configuration index $m$ then tells you how many of those will go through the north pole, and the rest will be judiciously spaced "parallels". If you insist on using single point charges, place alternating charges on each region, and the leading multipole contribution will be of order $l$.