## Find regions of stability of the Hill's equation
## Uses hill.m and hillstrob.m

function retval = hillstab(A, E) 
  empty_list_elements_ok = 1;
  retval = zeros(length(A), length(E));
  for k = 1:length(E)
    for j = 1:length(A)
      param = [A(j), E(k)];
      y1 = hillstrob([1, 0, param]);
      y2 = hillstrob([0, 1, param]);
      chpoly = poly([y1(2,1:2); y2(2,1:2)]);
      ## Condition of stability is | trace | <= 2
      retval(j, k) = chpoly(2);
    endfor
  endfor
endfunction

## ;;; Local Variables: ***
## ;;; mode: octave ***
## ;;; End: ***
