code: def model1(inps, lims, logger=None):
      
      a_arr, c_ = inps
      ilim, jlim = lims
  
      i_bs = sympy.symbols('i_lb i_ub', integer=True)
      i = sympy.Idx('i', i_bs)
  
      j_bs = sympy.symbols('j_lb j_ub', integer=True)
      j = sympy.Idx('j', j_bs)
  
      a = sympy.IndexedBase('a')
  
      c = sympy.Symbol('c', real=True)
  
      x = sympy.IndexedBase('x')
      y = sympy.IndexedBase('y')
  
      eqs = [
          sympy.Eq(x[i], (a[i]/3-1)**i+c),
          sympy.Eq(y[j], a[j]-j),
      ]
  
      ex_code = ExampleCode(eqs, (a[i], c), (i, j),
                            logger=logger, save_temp=True)
      x_, y_ = ex_code(inps, bounds=(ilim, jlim))
      x_ref = (a_arr/3-1)**np.arange(ilim[0], ilim[1]) + c_
      y_ref = a_arr[jlim[0]:jlim[1]] - np.arange(jlim[0], jlim[1])
      assert np.allclose(x_, x_ref)
      assert np.allclose(y_, y_ref)
docs: x[i] = (a[i]/3-1)**i + c
      y[j] = a[j] - j