Produce a matrix of b-splines for use internally within a Stan program

b_spline_function.stan

functions {
  // B function
  vector B(vector x, vector t, int i, int j_p1);
  
  vector B(vector x, vector t, int i, int k) {
    vector[rows(x)] out;
    vector[rows(x)] a_i1j1;
    vector[rows(x)] a_ij1;
    
    
    if(k==1) {
      for(n in 1:rows(x)) {
        if((t[i] <= x[n] && x[n] < t[i+1] )){ 
          out[n] = 1.0;
        } else {
          out[n] = 0.0;
        }
      }
    } else {
      if(t[i+k-1] == t[i]) {
        a_ij1 = rep_vector(0.0, rows(x));
      } else {
        a_ij1 = (x - t[i]) ./ (t[i+k-1] - t[i]);
      }
      
      if(t[i+k] == t[i+1]) {
        a_i1j1 = rep_vector(0.0, rows(x));
      } else {
        a_i1j1 = (t[i+k] - x) ./ (t[i+k] - t[i+1]);
      }
      
      out = a_ij1 .* B(x, t, i, k-1) + a_i1j1 .* B(x, t, i+1, k-1);
    }
    return(out);
  }
  
  // Create b-spline matrix
  matrix b_spline(vector x, vector interior_knots, int order) {
    real xmin;
    real xmax;
    vector[rows(interior_knots) + 2*order] t;
    matrix[rows(x), rows(interior_knots) + order] out;
    // fill out augmented knots vector
    xmin = min(x);
    xmax = max(x);
    
    t= append_row(append_row(rep_vector(xmin, order), 
                  interior_knots),
                  rep_vector(xmax, order));
    
    for(p in 1:(rows(interior_knots) + order)) {
      out[1:rows(x),p] = B(x, t, p, order-1);
    }
    for(i in 1:rows(out)) {
      if(x[i]==xmax) {
        out[i,cols(out)] = 1;
      }
    }
    return(out[1:rows(out),2:cols(out)]);
 }
}
data {
  
}
transformed data {

}
parameters {
  
}
transformed parameters {
  
}
model {
  
  
}