function [x,y] = stabdom(algor,typei)
%
%  Find the stability domain for an ODE solver
%
echo on
if typei == 1,
    alphai = 0;
else
    alphai = 180;
end
%
%  Loop over alpha
%
hvec = 0;
avec = 0;
count = 1;
terminate = -1;
alpha = alphai;
alphas = alphai;
type = typei;
dalpha = 5*typei;
hmn = 0;
hmx = 100;
swcount = 0;
echo off
while terminate < 0,
  hmax = hh(alpha,algor,hmn,hmx,type);
  if hmax > 0.0,
    alphas = alpha;
    avec = [ avec ; alpha ];
    hvec = [ hvec ; hmax ];
    count = count + 1;
    alpha = alphas + dalpha;
   else
    dalpha = dalpha/2;
    alpha = alphas + dalpha;
    if abs(dalpha) < 0.5,
      swcount = swcount + 1;
      type = -type;
      alpha = alphas;
      dalpha = 5*type;
    end,
    if swcount == 2,
      terminate = 1;
    end,
    if alpha > 180,
      terminate = 1;
    end,
    if alpha < 0,
      terminate = 1;
    end,
  end,
  hmn = 0;
  hmx = 1.5*hvec(count);
end
echo on
avec = [ avec ; 0 ];
hvec = [ hvec ; 0 ];
count = count + 1;
avec = [ avec(2:count) ];
hvec = [ hvec(2:count) ];
count = count - 1;
%
%  Convert to carthesian coordinates
%
radalpha = avec*pi/180;
x = - hvec .* cos(radalpha);
y = + hvec .* sin(radalpha);
%
%  Duplicate to get third quadrant also
%
xx = + x(count:-1:1);
yy = - y(count:-1:1);
x = [ x ; xx ];
y = [ y ; yy ];
%
return
