simplesys(nl) := flatten (makelist (ev ([
 -omega*epsilon[i] + h[i]*mu[i]*k = 0, 
 -h[i]*mu[i]*omega = -g*h[i]*sum(epsilon[i],i,0,nl-1)*k-h[i]*k*eta[i], 
 -h[i]*delta[i]*omega*%i = eta[i]-eta[i+1], 
  h[i]*mu[i]*k*%i+delta[i]-delta[i-1] = 0
                         ],
  eta[nl]: (eta[nl-2] - 6*eta[nl-1])/3, delta[-1]: 0, eta[-1]: eta[0]),
            i,nl-1,0,-1))$;

disdet(sys,nl) := eliminate(sys,
          delete(delta[0],flatten(makelist([epsilon[i],mu[i],delta[i],eta[i]],i,nl-1,0,-1))))[1]$;

dispersion(nl,sys) := map(rhs, solve(disdet(sys,nl), omega^2))[1]$;

omega2: makelist(dispersion(i,simplesys(i)),i,1,4);

kellersys(nl) := flatten(makelist (ev ([
 -omega*epsilon[i] + h[i]*mu[i]*k = 0, 
 -h[i]*mu[i]*omega = -g*h[i]*sum(epsilon[i],i,0,nl-1)*k-h[i]*k*(eta[i]+eta[i+1])/2, 
 -h[i]*(delta[i]+delta[i-1])/2*omega*%i = (eta[i]-eta[i+1]), 
  h[i]*mu[i]*k*%i+delta[i]-delta[i-1] = 0
                         ],
  delta[-1]: 0, eta[nl]: 0),
            i,nl-1,0,-1))$;

omega2_keller: makelist(dispersion(i,kellersys(i)),i,1,4);