I don't know where you got your formula from, but but I derived it this way:
Field inside the solenoid$=\mu_0ni \hat{z}$ (say)
Since the material is ferromagnetic, there is an induced, bound surface current $K\hat{\phi}$ (and $K=M$, where $M$ is magnetization). The magnetization is uniform, so bound current is zero,$$
J_b~=~\nabla \times \left(M\hat{x}\right)~=~0
\,.$$
From the Lorentz force, $F=i \vec{l}\times \vec{B}$:
$$
\begin{align}
\Rightarrow F &= A\vec{K}\times\vec{B} \\
&=AK\left(\mu_0ni\right) \left(\hat{\phi}\times\hat{z}\right) \\
&=AM(\mu_0ni)\hat{r} \\
&=A(\mu_0ni)(\chi_m*H)\hat{r} \\
&=A(\mu_0ni)(\chi_m*\frac{\mu_0ni}{\mu_0})\hat{r} \\
&=\chi_mA\mu_0 \left(ni \right)^2\hat{r}\,,
\end{align}
$$where $A$ is area.
Here $A$ is the area of the surface that $K$ is flowing on, i.e., the curved surface of the cylinder$=2\pi RL$ where $R$ is radius and $L$ the length of ferromagnet. The force is radially out of the surface of the core, stretching it out as if to fill the coil.
I don't know why the force should depend on the "gap" between the two.
(a) First we give an approximate force-based treatment for a solenoid whose length, $l$, is several times greater than its radius, $r$.
Consider the North end of the solenoid. There's a very simple argument to show that half the flux in the main part of the solenoid (that is half of $\pi r^2\ \mu_0 N I/l$) escapes through the sides of the solenoid before reaching the geometrical end of the solenoid. Therefore if $B_{rad}$ is the radial component of flux density at the surface of the solenoid at distance $z$ from its end,
$$\int_0^{z_0} B_{rad}\ 2\pi r dz =\tfrac 12 \frac{\pi r^2\ \mu_0 N I}l.$$
in which $z_0$ is the notional distance from the geometrical end at which there is no more radial escape of flux.
However in length $dz$ of solenoid there are $dz\ (N/l)$ turns, so a length $2\pi r\ dz\ (N/l)$ of wire. Therefore the total (axial) force on the North end part of the solenoid due to the radial field will be
$$\int_0^{z_0} B_{rad}\ \ I\ 2\pi r \tfrac {N}l dz\ \ =\ \tfrac {NI}l \int_0^{z_0} B_{rad}\ 2\pi r dz\ =\ \frac{\pi r^2\ \mu_0 N^2 I^2}{2l^2}.$$
This force acts towards the middle (lengthways) of the solenoid; a similar force acts on the South Pole end of the solenoid, so there is a compressive stress on the solenoid.
[It might be argued that, rather than $2\pi r$, the length of a turn is $\sqrt{(2\pi r)^2+(l/N)^2}$. So indeed it is, but the axial component, $(l/N)\hat z$, of the turn gives rise in the presence of $B_{rad}$ to tangential forces (that cancel over a complete turn anyway) rather than to an axial magnetic force!]
(b) The result just obtained agrees with what we get using $|F_{\text{axial}}|= -\frac {\partial U}{\partial l}$ in which $U$ is the magnetic energy stored in the solenoid. So, (ironically) ignoring end-effects and treating $B$ as uniform throughout the interior volume of the solenoid,
$$U=\frac 1{2\mu_0} B^2 \pi r^2 l=\frac 1{2\mu_0} \left( \frac {\mu_0 NI}l\right)^2 \pi r^2 l=\frac{\mu_0 N^2 I^2}{2l}\pi r^2.$$
$$\text{so}\ \ \ \ |F_{\text{axial}}|= -\frac {\partial U}{\partial l}=\frac{\mu_0 N^2 I^2}{2l^2}\pi r^2.$$
The irony is that in this method, (b), we have ignored what was the very basis of the force method, (a) – end effects! The reason it is permissible in this method is our assumption that $l>>r$. The end effects (escape of flux through the sides of the solenoid) are pretty much confined to end lengths $z_0$ of solenoid of the same order as the solenoid radius, $r$.
Best Answer
Yes, there are equations for the amount of force experienced by a ferromagnetic material in an external field - e.g. see Magnetic Force on a Ferromagenetic Material