Hi!

I defined some functions in a stan-mathematica interface to estimate parameters of a cosmological inhomogeneous model.

My problem is that i need some derivatives of those functions and to solve a differential equation to obtain the last parameters. For example i defined my most importan function (scale factor) in terms of other previous functions:

```
real[] R(real r, // Comovil Position
real t, // Comovil Time
real H0 , // Local Hubble Rate
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return (R0(r) / contrast(r, Omegaout, Omegain, r0, Deltar)) *
InverseSerie((2/3) * (contrast(r, Omegaout, Omegain, r0, Deltar) *
a(r, t, H0, Omegaout, Omegain, r0, Deltar))^(3/2));
}
```

And now i need to evaluate this function in some r(z) and t(z) values given by this system of equations that i could solve with mathematica, but i dont know how to implement in the stan code:

```
DtG[z_] :=
D[R[rG[z], tG[z]],
rG[z]]/((1 + z) D[D[R[rG[z], tG[z]], rG[z]], tG[z]])
DrG[z_] :=
Sqrt[1 - k[rG[z]]]/((1 + z) D[D[R[rG[z], tG[z]], rG[z]], tG[z]])
trsol = NDSolve[{tG'[z] == -DtG[z], rG'[z] == DrG[z], rG[0] == 0,
tG[0] == 1/H0}, {tG, rG}, {z, 0, 1.74}]
```

It gives me the values to perform the evaluation of the R (scale factor function) in the r(z), t (z) solutions and obtain R as function of z (and this is the function that i need to perform bayesian analisis)

Can i import the function that i have obtained in mathematica via NDSolve? There is a way to implement this in the â€śStanCodeâ€ť interface for mathematica? I read some information about autodiff useful, but i dont know how to apply in my case.

I would very grateful if someone can give me an orientation cause im a little lost and i like stan as bayesian MC estimator!!

Thanks!