Mathematically speaking, it is true that Maxwell's equations by themselves aren't the whole story; they're a set of PDEs for which one needs to specify boundary conditions separately if one wants to solve them. Boundary conditions can be well-motivated from a physical perspective in a given scenario, but they do not follow from the equations themselves.
As for Poisson's equation versus Coulomb's Law, no requirement of spherical symmetry is necessary. Start with Poisson's equation, and set the charge density to be that of a point charge, namely
\begin{align}
\rho(\mathbf x) = q\delta(\mathbf x - \mathbf x_0).
\end{align}
Next, use the vanishing curl of the electric field (which works in the absence of explicitly time-varying magnetic fields) to write $\mathbf E = -\nabla \Phi$ and plug this into Poisson's equation to obtain
\begin{align}
\nabla^2\Phi = -\frac{q}{\epsilon_0}\delta(\mathbf x-\mathbf x_0)
\end{align}
In other words, we want to determine $\Phi$ that is a Green's function for the Laplace equation. The general solution is
\begin{align}
\Phi = \frac{1}{4\pi\epsilon_0}\frac{q}{|\mathbf x - \mathbf x_0|} + F(\mathbf x)
\end{align}
where $F$ is a harmonic function, namely one which satisfies Laplace's equation. If we then invoke the boundary condition that the potential vanishes at infinity, then this forces $F$ to be identically zero, and we obtain
\begin{align}
\Phi = \frac{1}{4\pi\epsilon_0}\frac{q}{|\mathbf x - \mathbf x_0|}
\end{align}
which, upon taking the gradient, yields he electric field of a point charge which is essentially the content of Coulomb's Law.
The auxiliary magnetic field is not zero in your example. The reason is that $\mathbf H$ must also obey continuity conditions at the boundary and these act as source terms. These source terms are the surface equivalents of the 'magnetization charge' $\nabla\cdot\mathbf M$, and they're present at surfaces where the magnetization has a discontinuity $\Delta\mathbf M$ with a nonzero component across the surface normal $\hat{\mathbf n}$.
As always, the continuity conditions at the boundary the limiting form of Maxwell's equations at the boundary. In this case they come from the magnetic Gauss law and the magnetostatic Ampère law, which for the $\mathbf B$ field read
$$
\oint_S\mathbf B\cdot\mathrm d\mathbf S=0
\quad\text{and}\quad
\oint_C \mathbf B\cdot \mathrm d\mathbf l=\mu_0I_\mathrm{enc},
\tag 1
$$
and in their boundary form
$$
\hat{\mathbf n}\cdot\Delta \mathbf B=0
\quad\text{and}\quad
\hat{\mathbf n}\times\Delta \mathbf B=\mu_0\mathbf K.
\tag 2
$$
This means that the normal component of $\mathbf B$ should be continuous across any surface. It should then be clear that the same cannot be true for the auxiliary $\mathbf H$ field, because $\mathbf H=\tfrac{1}{\mu_0}\mathbf B-\mathbf M$ is the sum of a continuous $\mathbf B$ field and a discontinuous magnetization $\mathbf M$.
Moreover, this discontinuity can be directly quantified from the boundary-form Maxwell equation by simply plugging in the definition of $\mathbf H$, to get
$$
\hat{\mathbf n}\cdot\Delta \mathbf H=-\frac{1}{\mu_0}\hat{\mathbf n}\cdot\Delta \mathbf M.
\tag 3
$$
This acts as a source term for $\mathbf H$ and it's what makes it nonzero throughout.
On the other hand, the specific case you mention, in which the entire space is filled with the magnetized material, is trickier to deal with. In general, infinite bulk sources (as opposed to infinite planes or lines) are iffy propositions in electromagnetism, and they don't always give physical results. (For an extreme example, a constant charge density throughout all space, there's simply no acceptable resolution.)
In such cases what one needs to do is to take a finite volume and then take the limit as it becomes bigger and bigger. One then needs to check that the answer is (i) independent of the shape of the volume, and (ii) independent of the relative position of the probe point with respect to the centre of expansion; if the answer is yes on both counts then that is your answer.
For the specific case when all of space is filled with a magnetized material, this seems to be doable. The easiest way to approach this is to use slabs of infinite area and finite width along the magnetization direction, and then take the width to infinity. Here in every finite case $\mathbf H=\mathbf M$ is easily verified from the boundary conditions, so $\mathbf B=2\mu_0\mathbf M$ is independent of the centre of expansion and therefore the correct answer.
As an independent check, it is relatively easy to verify that a constantly magnetized sphere of material produces a uniform magnetic field inside it, so you could take the limit as the radius goes to infinity and you will (presumably) get the same answer.
To round things up, I'll refer you to Jackson's Classical Electrodynamics §5.8 (pp. 191-194 on the 3rd edition, Wiley 1999), which provides a more thorough treatment. In particular, to get a complete formulation in terms of $\mathbf H$ you also need to provide a condition for its tangential component. To do this, you simply rephrase Ampère's law as
$$
\oint_C \mathbf H\cdot \mathrm d\mathbf l=I_\mathrm{free},
\tag 4
$$
which then boils down to
$$
\hat{\mathbf n}\times\Delta \mathbf H=\mathbf K_\mathrm{free}.
\tag 5
$$
In your case there are no free currents, so the tangential component of $\mathbf H$ is continuous everywhere (i.e. it behaves exactly like a static electric field!), but there is no escaping the discontinuity in the normal component (which is again analogous to that of an electric field).
Finally, note that Jackson includes some additional simplifications which can be brought to bear in linear media, in which case $\mathbf B=\mu\mathbf H$ and you can rephrase $(3)$ as $\hat{\mathbf n}\cdot \mathbf H_1=\frac{\mu_2}{\mu_1}\hat{\mathbf n}\cdot \mathbf H_1$ (assuming isotropic media, though that can be relaxed easily). However, this does not apply in your case since you have a fixed magnetization, so your medium is not linear.
Best Answer
Any static field with zero divergence, i.e. which obeys the magnetic Gauss law $\nabla\cdot\mathbf B=0$, is a valid magnetic field. The curl of the field can be anything: if it is nonzero then it requires a current density $\mathbf J$ to sustain it, given by Ampère's law $$ \nabla\times\mathbf B = \mu_0\mathbf J. $$ This curl can in principle be any suitable vector field, though because it is a curl it is required to have zero divergence, i.e. $\nabla\cdot(\nabla \times \mathbf B )= \mu_0\nabla \cdot\mathbf J=0$, which needs to happen anyway in a static situation, because of conservation of charge.
If there are no currents, i.e. in vacuum, then yes, the magnetic field will have zero curl. Most of the usual examples of magnetic fields fall into this category, and it is plenty possible for a magnetic field to have zero divergence and zero curl (want a simple example? try a constant field).
It is important to note, though, that typically you want your magnetic field to have sources somewhere, and this means that the in-vacuum condition $\nabla\times\mathbf B=0$ will only hold for some restricted region in space. In this spirit, it is possible to have a vector field with no divergence and no curl defined over all of space, but this requires its magnetic energy $U=\frac{\mu_0}{2}\int|\mathbf B|^2\mathrm d\mathbf r$ to be infinite, which is as unphysical as a nonzero static magnetic field with no sources.
If you want a time-dependent field, though, the problem changes a good bit. As far as the magnetic Maxwell equations are concerned, $$ \nabla\cdot\mathbf B=0 \text{ and } \nabla\times\mathbf B = \mu_0\varepsilon_0 \frac{\partial \mathbf E}{\partial t}+\mu_0\mathbf J, $$ any vector field with no divergence can in principle be interpreted as a magnetic field. However, this requires us to find a set of sources, $\mathbf J$ and $\rho$, and more importantly an electric field $\mathbf E(\mathbf r,t)$, to go with it, or the whole thing is a bit moot. This means, therefore, that we need to turn Maxwell's equations on their head a little bit, and now they become \begin{align} \mu_0\varepsilon_0 \frac{\partial \mathbf E}{\partial t}+\mu_0\mathbf J & = \nabla\times\mathbf B \\ \nabla\times\mathbf E & = - \frac{\partial \mathbf B}{\partial t} \end{align} as a problem to be solved for $\mathbf E$ and $\mathbf J$, which is generally always solvable.