If model evaluates at initial conditions without autodiff but fails with autodiff then the error messages are not great

Summary:

For instance, with an ODE model @jtimonen found output like this:

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: ode_bdf_tol:  Failed to integrate to next output time (1.542) in less than max_num_steps steps (in 'lv.stan', line 34, column 2 to line 35, column 48)
Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: ode_bdf_tol:  Failed to integrate to next output time (1.542) in less than max_num_steps steps (in 'lv.stan', line 34, column 2 to line 35, column 48)
Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: ode_bdf_tol:  Failed to integrate to next output time (1.542) in less than max_num_steps steps (in 'lv.stan', line 34, column 2 to line 35, column 48)
Exception: ode_bdf_tol:  Failed to integrate to next output time (1.542) in less than max_num_steps steps (in 'lv.stan', line 34, column 2 to line 35, column 48)
Exception: ode_bdf_tol:  Failed to integrate to next output time (1.542) in less than max_num_steps steps (in 'lv.stan', line 34, column 2 to line 35, column 48)

We would expect an error message like:

Rejecting initial value:
  Error evaluating the log probability with gradients at the initial value.
Exception: ode_bdf_tol:  Failed to integrate to next output time (1.542) in less than max_num_steps steps (in 'lv.stan', line 34, column 2 to line 35, column 48)

Initialization between (-2, 2) failed after 100 attempts. 
 Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
Initialization failed.

Here is a model that exhibits this behavior at certain seeds (model below):

functions {
  // Lotka-Volterra system                                                                                                                                                                    
  vector derivative_fun(real t, vector y, data int[] a0, vector theta) {
    vector[2] dydt;
    dydt[1] = theta[1]*y[1] - y[1]*y[2];
    dydt[2] = y[1]*y[2] - theta[2]*y[2];
    return dydt;
  }
}

data {
  int<lower=1> N;
  real t_eval[N]; // must be increasing                                                                                                                                                       
  vector[2] y_data[N];
  vector[2] y0;
  real t0;
  real<lower=0> ABS_TOL;
  real<lower=0> REL_TOL;
  int<lower=1> MAX_NUM_STEPS;
}

transformed data {
  int a0[0];
}

parameters {
  vector<lower=0>[2] theta;
  real<lower=0> sigma;
}

model {
  theta ~ normal(1, 0.3);
  sigma ~ normal(0, 2.0);
  vector[2] y_hat[N] = ode_bdf_tol(derivative_fun, y0, t0, t_eval,
    REL_TOL, ABS_TOL, MAX_NUM_STEPS, a0, theta);
  for(n in 1:N){
    target += normal_lpdf(y_data[n] | y_hat[n], sigma);
  }
}

Json data:

{
  "N": 20,
  "t_eval": [1.542, 3.314, 3.582, 3.942, 4.618, 5.974, 6.338, 6.546, 6.698, 7.106, 7.53, 7.602, 7.758, 8.334, 8.45, 9.018, 9.17, 10.886, 11.366, 11.574],
  "y0": [1, 2],
  "t0": 0,
  "y_data": [
    [0.344797088393873, 1.6930242060351],
    [0.635418289383446, 0.0906638098070269],
    [0.519516090322978, 0.821024324067281],
    [1.07889503979525, 0.477741038303985],
    [0.773776915048449, 0.255210793923543],
    [1.72581484166438, 1.5527620189067],
    [0.851008213891165, 1.80589311943048],
    [1.62769164960979, 1.57197163973671],
    [1.06907814440105, 1.77751680937944],
    [0.559461058239121, 1.77778520957167],
    [0.319224257399646, 2.09229342823705],
    [0.656623356899921, 2.35616034285975],
    [0.266800987984575, 1.59088141802012],
    [-0.0843099951671533, 1.29952132205262],
    [0.510139049085921, 1.16569483217513],
    [0.739499156159327, 1.5954691564362],
    [0.300537914354303, 1.33024642117391],
    [0.453591698548004, 0.627159819286871],
    [1.13100675948511, 0.372307762827462],
    [1.63136411638082, 1.04814114632482]
  ],
  "REL_TOL": 1e-06,
  "ABS_TOL": 1e-06,
  "MAX_NUM_STEPS": 30
}

Current Version:

v2.26.1

Start at the ode_bdf_tol entry point and reproduce the Lotka-Volterra model from the issue using the supplied JSON data and v2.26.1 settings. Compare initialization failures without autodiff against failures with gradients, including the messages shown for lv.stan. Done means gradient-aware wording is reported and the expected initialization-failure guidance appears without repeated misleading messages.