ODE question

Base Code 10.26 in the book of Bayesian Models for Astrophysical Data Using R, JAGS, Python, and Stan for jla data :

functions {

real[] Ez(real z,
real[] H,
real[] params,
real[] x_r,
int[] x_i) {
real dEdz[1];
dEdz[1] = 1.0/sqrt(params[1]*(1+z)^3
return dEdz;
data {
int<lower=1> nobs; // number of data points
real E0[1]; // integral(1/H) at z=0
real z0; // initial redshift, 0
real c; // speed of light
real H0; // Hubble parameter
vector[nobs] obs_mag; // observed magnitude at B max
real x1[nobs]; // stretch
real color[nobs]; // color
real redshift[nobs]; // redshift
real hmass[nobs]; // host mass
transformed data {
real x_r[0]; // required by ODE (empty)
int x_i[0];
real<lower=0, upper=1> om; // dark matter energy density
real alpha; // stretch coefficient
real beta; // color coefficient
real Mint; // intrinsic magnitude
real deltaM; // shift due to host galaxy mass
real<lower=0> sigint; // magnitude dispersion
real<lower=-2, upper=0> w; // dark matter equation of state
transformed parameters{
real DC[nobs,1]; // co-moving distance
real pars[2]; // ODE input = (om, w)
vector[nobs] mag; // apparent magnitude
real dl[nobs]; // luminosity distance
real DH; // Hubble distance = c/H0
pars[1] = om;
pars[2] = w;
DH = (c/H0);
# Integral of 1/E(z)
DC = integrate_ode_rk45(Ez, E0, z0, redshift, pars, x_r, x_i);
for (i in 1:nobs) {
dl[i] = DH * (1 + redshift[i]) * DC[i, 1];
if (hmass[i] < 10) mag[i] = 25 + 5 * log10(dl[i])
+ Mint - alpha * x1[i] + beta
* color[i];
mag[i] = 25 + 5 * log10(dl[i])
+ Mint + deltaM - alpha * x1[i] + beta * color[i];
model {
# Priors and likelihood
sigint ~ gamma(0.001, 0.001);
Mint ~ normal(-20, 5.);
beta ~ normal(0, 10);
alpha ~ normal(0, 1);
deltaM ~ normal(0, 1);
obs_mag ~ normal(mag, sigint);

We can use this code when we have Analytical solutions and put these solutions in Function block
such as:\Omega_{D}(z)= (1-params[1])(1+z)^{(3(1+params[2]))} or params[1](1+z)^{3}. My question is
what if this parameter \Omega_{D} is solution of ode such as
\Omega_{D}'=\frac{d\Omega_{D}}{d z}=-3(\delta-1)\Omega_{D}\frac{1-\Omega_{D}}{1-(2-\delta)\Omega_{D}} that \delta
and initial value for \Omega_{D} are our parameters .
For a better understanding of problem and what do I want Solve problem with python.

def ode_func(x, t,a):
    return -(3*pow(1+t,-1)*(a-1))*x*((1-x)/(1-2*x+a*x))

sol=lambda a,inits: odeint(ode_func,inits,t,args=(a,)) 

def model(a,inits):
    y = np.ravel(sol(a,inits))
    return interp1d(t, y, kind='cubic')

def comoving_integrand(z, omm,a):
    f2= model(a,inits)
    Ez = 1.0/np.sqrt((omm*(1.0+z)**3 + f2(z)))
    return Ez
def luminosityDistance(z, omm,a):
    integral = integrate.quad(comoving_integrand, 0, z, args=(omm,a))[0]
    return 300000*(1+z)*integral 

My goal is achieving luminosityDistance that is similar to DC = integrate_ode_rk45(Ez, E0, z0, redshift, pars, x_r, x_i);
in transformed parameters block of Stan code that I mentioned above .
Any suggestions would be greatly appreciated.

I don’t see a problem with that…

are you having trouble specifying the ODE function? the model? having inference issues?

I don’t know how define function that first solve ODE for \Omega_{D} and then use this solution into the dEdz function block as I wrote the python script, first solve ODE for \Omega_{D} then interpolation over interval give me f2(z) that used as integrate function comoving_integrand(z, omm,a).