ODE with random effects running slow

Modeling
1

Hi all,

After running my ODE with a random effect for the initial population, I get a computational time of 0.455 seconds but it still goes extremely slow. I am running this with stan, but I have also run this with cmdstan which did work without any warning.
My initial parameters look like this:

initebl6 = function(){
list(ab = -0.4, logAB0 = log(36), sigmalogAB0 = 0.5,
rAB0 = rep(0, participants),
sigma = 0.5)
}

Also, my data consists of 700 individuals with 3 samples per individual, so my times go per individual from 0,57, 78 with some variability within these.

Do you have any suggestions that I should try on? Why is this model running this slow when the gradient evaluation is not technically large?

 functions {
  real[] edecay(real t, real[] y, real[] theta,  real[] x_r, int[] x_i) {
    real dydt[1];
    dydt[1] = -theta[1]*y[1];
    
    return dydt;
  }
}
data {
  int<lower=1> N;
  int<lower=1> nsubjects;
  int<lower=1> subj_ids[N];
  real <lower=0> y[N];
  int <lower=1> ntimes;
  real ts[ntimes];
  real t0;
}
transformed data {
  real x_r[0];
  int x_i[0];
}
parameters {
  real ab;
  //real <lower = 0> sigmalogab;
  real <lower=0> logAB0;
  real  sigmalogAB0;
  vector[nsubjects] rAB0;
  real <lower =0> sigma;
}
model{
  //prior distributions
  real new_mu[N, 1];
  real mu[N, 1];
  real inits[1];
  real theta[1];
  //real AB0 = exp(logAB0);
  real sigmaAB0 = exp(sigmalogAB0);
  real eval_time[ntimes];
  theta[1] = exp(ab);
  //prior distributions
  ab ~ normal(0, 1);
  //sigmalogmuab ~ exponential(1.0);
  logAB0 ~ lognormal(log(36), 1);
  sigmalogAB0 ~ normal(0, 0.5);
  sigma ~ normal(0, 1);
  //rho_rhobit ~ uniform(-10, 10);
  for (j in 1:nsubjects){
    rAB0[j] ~ normal(-sigmaAB0/2, sigmaAB0);
  }
  for( i in 1:N){
  inits[1] = logAB0*exp(rAB0[subj_ids[i]]);
  eval_time[i] = ts[i];
  if (eval_time[i] == 0){
    new_mu[i,1] = log(inits[1]) - sigma/2;}
  else{
  mu[i:i] = integrate_ode_rk45(edecay, inits, t0, eval_time[i:i], theta, x_r, x_i); 
  new_mu[i,1] = log(mu[i, 1]) - sigma/2;
  }
  y[i] ~ lognormal(new_mu[i, 1], sigma);
  }
}

Thanks in advance.