cciopi.py wrote: "$\rm \mathcal H=a'(\eta)/a(\eta)$ where $\rm a'(\eta)$ is the derivative of the scale factor with respect to the conformal time."
This expression for $\mathcal H$ is just the reciprocal of the comoving Hubble radius $\rm \mathcal H = H a$, so that is trivial since if you have $\rm a$ you also know $\rm H$, so you just need to multiply them. More interesting is $\rm H$ in terms of $\eta$, so I'll answer that as well although you didn't ask for it.
cciopi.py wrote: "Specifically, I'd like to know how to compute $\rm \mathcal H$
from the cosmological model, such as the ΛCDM model, for a given redshift $\rm z$."
I suppose you already know the relation $\rm a=1/(z+1)$. The particle horizon $\rm r_p$ in proper distances is the integral
$$\rm \rm r_p = a \int_0^a \frac{c}{\bar a^2 \ H(\bar a)} d \bar a = a(t)\int_0^t \frac{c}{a(\bar t)} d \bar t$$
and $\rm \eta = R_p/c$ where $\rm R_p=r_p/a$ is the particle horizon in comoving distances.
so for example if you take a radiation dominated universe where $\rm \Omega_R=1$
$$\rm H(a)=\frac{H_0}{a^2} \ \to \ a(t)=\sqrt{2 \ H_0 \ t} \ \to \ H(t)=\frac{1}{2 \ t}$$
the relation of $\rm r_p$ in terms of $\rm a$ or $\rm t$ is
$$\rm r_p = \frac{c \ a^2}{H_0} = 2 \ c \ t \ \to \ t=\frac{r_p}{2 \ c}$$
which gives the conformal time $\rm \eta$ in terms of $\rm a$ or $\rm t$
$$\rm \eta = \frac{a}{H_0} = \sqrt{\frac{2 \ t}{H_0}} \ \to \ t=\frac{H_0 \ \eta^2}{2}$$
or in a matter dominated universe where $\rm \Omega_M=1$
$$\rm H(a)=\frac{H_0}{\sqrt{a^3}} \ \to \ a(t)=\left(\frac{3 \ H_0 \ t}{2}\right)^{2/3} \ \to \ H(t)=\frac{2}{3 \ t}$$
the relation in terms of $\rm a$ or $\rm t$ is
$$\rm r_p = \frac{2 \ c \ \sqrt{a^3}}{H_0} = 3 \ c \ t \ \to \
t=\frac{r_p}{3 \ c}$$
which gives the conformal time $\rm \eta$ in terms of $\rm a$ or $\rm t$
$$\rm \eta = \frac{2 \sqrt{a}}{H_0} = \sqrt[3]{\frac{12 \ t}{H_0^2}} \ \to \
t=\frac{H_0^2 \ \eta^3}{12}$$
or in a dark energy dominated universe where $\rm \Omega_{\Lambda}=1$
$$\rm H(a)=H_0 \ \to \ a(t)=e^{ \ H_0 \ t} \ \to \ H(t)=H_0$$
the relation in terms of $\rm a$ or $\rm t$ is
$$\rm r_p = \frac{c \ (a-1)}{H_0} = \frac{c \ (e^{ \ H_0 \ t}-1)}{H_0} \ \to \ t=\frac{\ln(1+H_0 \ r_p/c)}{H_0}$$
which gives the conformal time $\rm \eta$ in terms of $\rm a$ or $\rm t$
$$\rm \eta = \frac{1-1/a}{\ H_0} = \frac{1-e^{-H \ t}}{H} \ \to \ t=\frac{\ln(\frac{1}{1-\eta \ H_0})}{H_0}$$
so plug the $\rm a$ or $\rm t$ of $\rm \eta$ into the equation for $\rm H(a)$, $\rm H(t)$ or whatever you need as a function of $\rm \eta$. For radiation domination this gives $\rm H=1/(H_0 \ \eta^2)$ and $\rm \mathcal H=1/\eta$, for matter domination $\rm H=8/(H_0^2 \ \eta^3)$ and $\rm \mathcal H=2/\eta$ etc.
In the last example with dark energy domination you have to integrate $\rm \int_1^a$ instead of $\rm \int_0^a$ since $\rm a$ is never zero in such an universe, so we usually start with $\text{a=1}$ at $\text{t=0}$ in that case. For $\rm H$ this makes no difference since it's constant in that case anyway, but you might need it for the other functions. $\rm \mathcal H=H_0/(1-H_0 \ \eta)$ in that case though if you really need it.
I don't know about your program, but for my code see here. For different or mixed component universes see here for the corresponding $\rm H(t)$ or $\rm a(t)$.
In our universe with a combination of $\rm \Omega_R$, $\rm \Omega_M$ & $\rm \Omega_{\Lambda}$ there are no analytical solutions and you have to do it numerically.
If you neglect $\rm \Omega_R$ and only take $\rm \Omega_M$ & $\rm \Omega_{\Lambda}$ there is an analytical solution though, but it is not really elegant, so I won't paste it here, but you can find it in the last link above.
