function plotS(R,plotOS,axisbox)
%
% plot S, the region of absolute stability for the rational function R(z).
%
% If plotOS=1 then also plot the Order Star (complement of the 
% region of relative stability) 
%
% Typically R(z) is an approximation to exp(z) obtained by applying a 
% finite difference method to y' = lambda y, where z = dt*lambda.
%
% axisbox = axis region 
%           (if omitted, default values are used that may be a poor choice)
%
% See also makeplotS.m, from which this routine is called for a variety of
% specific R(z) functions.
%
% From  http://www.amath.washington.edu/~rjl/fdmbook/chapter7  (2007)

if nargin==1     % if only one input argument R:
   plotOS = 0;   % defaults to plotting only S.
   end

if nargin==3
    xa = axisbox(1); xb = axisbox(2); ya = axisbox(3); yb = axisbox(4);
  else
    xa = -10; xb = 10; ya = -10; yb = 10;  % default values
  end

% compute ratio R(z) over a fine grid on the region [xa,xb] x [ya,yb]

nptsx = 501;
nptsy = 501;
x = linspace(xa,xb,nptsx);
y = linspace(ya,yb,nptsy);
[X,Y] = meshgrid(x,y);
Z = X + 1i*Y;
Rval = feval(R,Z);

%----------------------------------------------------------------------
% Plot S, region of absolute stability:

Rabs = abs(Rval);

% plot Rabs with a color map that is shaded in S (where Rabs < 1)
% and which in complement (where Rabs > 1).

figure(1)
clf
surf(x,y,0*Rabs,Rabs,'EdgeColor','none')
view(2)
colormap([.7 .7 1;1 1 1])
%colormap([.7 .7 .7;1 1 1])   % grayscale
caxis([.99 1.01])
daspect([1 1 1])
hold on
contour(x,y,Rabs,[1 1],'k')
box on
hold off
title('Region of absolute stability','FontSize',15)
% plot axes:
hold on
plot([xa xb],[0 0],'k')
plot([0 0],[ya yb],'k')
axis([xa xb ya yb])
set(gca,'FontSize',15)
hold off

if ~plotOS, return, end

%----------------------------------------------------------------------
% Plot order star, complement of region of relative stability:

Rrel = abs(Rval ./ exp(Z));

% plot Rrel with a color map that is shaded in order star (where Rrel > 1)
% and white in complement (where Rrel < 1).

figure(2)
clf
surf(x,y,0*Rrel,Rrel,'EdgeColor','none')
view(2)
colormap([1 1 1; .7 .7 1])
%colormap([1 1 1; .7 .7 .7])   % grayscale
caxis([.99 1.01])
daspect([1 1 1])
hold on
contour(x,y,Rrel,[1 1],'k')
box on
hold off
title('Order Star','FontSize',15)
% plot axes:
hold on
plot([xa xb],[0 0],'k')
plot([0 0],[ya yb],'k')
axis([xa xb ya yb])
set(gca,'FontSize',15)
hold off