function h = stabh(D,s)

%  stabh - Determines largest stable time step for an s order Runge-Kutta
%  method (s<5) on y' = D*y, with coeffient matrix D. 


  if nargin < 2
    error('Not enough input arguments.');
  end

  if nargin > 2
    error('Too many input arguments.');
  end
  
  h = 1;
  lam = eig(D); % calculate eigenvalues for testing stability region
  i = 0;    
  while any(Rh(h*lam,s) >= 1.0),
    h = h/2;  % step size a power of 2 for less roundoff on computers
    i = i+1;  % count order of h
    if i >= 30 % if h gets too small, exit
      error('No suitable stepsize h found (h < 1e-9 required).');
    end
  end
  
function r = Rh(hl,s)

  r = hl + 1; 
  
  fac = 1;
  for i=2:s,
    fac = fac*i;  %i factorial
    r = r + hl.^i/fac;
  end
  r = abs(r);

  
