Hello,

I ran into a puzzling error when testing a multivariate normal model and was hoping you could help me understand what the issue is. Basically, the same covariance function written in two different forms, one will run and the other will not, even though they are exactly the same function. What am I missing? I paste the code for a simplified version of the model that reproduces this issue. Data are attached (n = 2). I am using cmdstanpy (stan version 2.23).

```
functions{
matrix matern52_v1(matrix qx,matrix qy, real s, real l){
matrix[rows(qx),cols(qx)] r = sqrt(qx + qy)/l;
return square(s)*(1 + sqrt(5)*r + 5/3*(r .* r)) .* exp(-sqrt(5)*r);
}
matrix matern52_v2(matrix qx,matrix qy, real s, real l){
matrix[rows(qx),cols(qx)] r = sqrt(qx/square(l) + qy/square(l));
return square(s)*(1 + sqrt(5)*r + 5/3*(r .* r)) .* exp(-sqrt(5)*r);
}
}
data {
int<lower=0> n;
vector[n] lon;
vector[n] lat;
vector[n] y;
}
transformed data{
matrix[n,n] qx;
matrix[n,n] qy;
for(i in 1:n){
for(j in 1:n){
qx[i,j] = (lon[i] - lon[j]) * (lon[i] - lon[j])/100;
qy[i,j] = (lat[i] - lat[j]) * (lat[i] - lat[j])/100;
}
}
}
parameters {
real<lower=0> l;
real<lower=0> s;
}
model {
matrix[n,n] q;
matrix[n,n] L;
l ~ std_normal();
s ~ std_normal();
q = matern52_v1(qx,qy,s,l);
//q = matern52_v2(qx,qy,s,l);
L = cholesky_decompose(q);
y ~ multi_normal_cholesky(rep_vector(0,n),L);
}
```

Note that the functions ‘matern52_v1’ and ‘matern52_v2’ are exactly the same function just written differently. However, when I use matern52_v1 the code will run fine whereas when I use matern52_v2 the code will throw the following error:

```
RuntimeError: Error during sampling.
chain 1 returned error code 70
```

The *stdout.txt file contains the following message:

```
Rejecting initial value:
Gradient evaluated at the initial value is not finite.
Stan can't start sampling from this initial value.
```