function hfig = draw_stabfunc(f,ReMin,ReMax,ImMin,ImMax,N)
% stability_region Plot the stability region of a function f
%
%   INPUT:      f       Function of two variables a and b, the real and
%                       imaginary parts of h*lmax. Returns the compression
%                       factor for any combination of a and b, must be
%                       vectorized.
%               ReMin   Mininum value for the real part of h*lmax
%               ReMax   Maximum value for the real part of h*lmax
%               ImMin   Mininum value for the imaginary part of h*lmax
%               ImMax   Maximum value for the imaginary part of h*lmax
%               N       Number of plot points
%
%   OUTPUT:     hfig    Figure where the plot resides
%
%   EXAMPLE:    
%       The Euler Forward stability region is given by abs(1+h*lmax) < 1.
%       Therefore the function will be f = @(a,b)sqrt((1+a).^2 + b.^2);
%       The plot can now be generated by using
%       hfig = draw_stabfunc(f,-2.5,0.5,-1.5,1.5,1000);

ha = linspace(ReMin,ReMax,N);
hb = linspace(ImMin,ImMax,N);

[X,Y] = meshgrid(ha,hb);

Z = f(X,Y);

Z(Z>=1) = nan;

hfig = figure;
contourf(X,Y,Z,0)
colormap([0.5000    1.0000    0.5000])

hold on
line([0,0],[min(hb),max(hb)],'linestyle','--','color','k','linewidth',2)
line([min(ha),max(ha)],[0,0],'linestyle','--','color','k','linewidth',2)
hold off

xlabel('h*Re(\lambda)','fontsize',12)
ylabel('h*Im(\lambda)','fontsize',12)
pbaspect([1, 1, 1])