
# function VL(var)
#  local res = {}
#  for i = 0, N-1 do
#    res[i+1] = var..'_'..i
#  end
#  return table.concat(res,',')
# end

# function AL(var)
#  local res = {}
#  for i = 0, N-1 do
#    res[i+1] = var..'['..i..']'
#  end
#  return table.concat(res,',')
# end

# function VLI(var, ord)
#  local res = {}
#  for i = 0, N-1 do
#    res[i+1] = string.format('%s%i_%i', var, ord, i)
#  end
#  return table.concat(res,',')
# end

# function KCONV(var, ord, i)
#    local sm = {}
#    for j = 1, ord do
#       local bc = var[j]
#       if tonumber(bc) ~= 0 then
#          sm[#sm+1] = string.format('(%s)*k%i_%i', bc, j, i)
#       end
#    end
#    return table.concat(sm, ' + ')
# end

local ffi  = require 'ffi'

local function ode_new()
   local n = $(N)
   return {t = 0, h = 1, dim = n, y = matrix.new(n, 1), dydt = matrix.new(n, 1)}
end

local function ode_init(s, t0, h0, f, $(VL'y'))
   $(AL's.y.data') = $(VL'y')
   $(AL's.dydt.data') = f(t0, $(VL'y'))
   s.t, s.h, s.f = t0, h0, f
end

local function ode_evolve(s, t1, tsmp)
   local step, t0, y = s.step, s.t, s.y

   local function it(s, t)
      t = t and t + tsmp or t0
      if t <= t1 then
	 while s.t < t do
	    step(s, t)
	 end
	 return t, $(AL'y.data')
      end
   end

   return it, s
end

local function hadjust(rmax, h)
   local S = 0.9
   if rmax > 1.1 then
      local r = S / rmax^(1/$(order))
      r = max(0.2, r)
      return r * h, -1
   elseif rmax < 0.5 then
      local r = S / rmax^(1/($(order)+1))
      r = max(1, min(r, 5))
      return r * h, 1
   end
   return h, 0
end
