I now have code that works:
functions {
// Example integrand
real[] integrand_ode(real y, real[] f, real[] theta, real[] x, int[] x_i) {
real df_dx[1]; // Do not change
df_dx[1] = exp(-y^2 - x[1]^2); // Function of x and y to be integrated
return(df_dx);
}
//User-defined function of x after y integrated out
real[] first_integral(real x, real[] f, real[] theta, real[] x_r, int[] x_i) {
return(integrate_ode_rk45(integrand_ode,
rep_array(0.0, 1), // Do not change
0.0, // ymin
rep_array(1.0, 1), //(ymax, 1)
rep_array(0.0, 0), // Do not change
rep_array(x, 1), // Do not change
x_i)[1,]);// Do not change
}
// Final integration wrt x
real second_integral(real xmin, real xmax ) {
int x_i[0];
return(integrate_ode_rk45( first_integral,
rep_array(0.0, 1),// Do not change
xmin, //xmin
rep_array(xmax, 1), // (xmax, 1)
rep_array(0.0, 0),// Do not change
rep_array(0.0, 0) , // Do not change
x_i)[1,1]);// Do not change
}
}
data {
}
model {}