function varargout = stability_gui(varargin)
% STABILITY_GUI MATLAB code for stability_gui.fig
%      STABILITY_GUI, by itself, creates a new STABILITY_GUI or raises the existing
%      singleton*.
%
%      H = STABILITY_GUI returns the handle to a new STABILITY_GUI or the handle to
%      the existing singleton*.
%
%      STABILITY_GUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in STABILITY_GUI.M with the given input arguments.
%
%      STABILITY_GUI('Property','Value',...) creates a new STABILITY_GUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before stability_gui_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to stability_gui_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help stability_gui

% Definitions of various Linear Multistep and One Step Methods

global lmm;

% ***LMM methods***
% r = rho
% s = sigma

% Backward Differentiation Formulae (BDF)

BDF1.r=[1 -1]; % AKA
BDF1.s=[1 0];  % backward Euler

BDF2.r=[1 -4/3 1/3];
BDF2.s=[2/3 0 0];

BDF3.r=[1 -18/11 9/11 -2/11];
BDF3.s=[6/11 0 0 0];

BDF4.r=[1 -48/25 36/25 -16/25 3/25];
BDF4.s=[12/25 0 0 0 0];

BDF5.r=[1 -300/137 300/137 -200/137 75/137 -12/137];
BDF5.s=[60/137 0 0 0 0 0];

BDF6.r=[1 -360/147 450/147 -400/147 225/147 -72/147 10/147];
BDF6.s=[60/147 0 0 0 0 0 0];

% Adams-Bashforth

AB1.r = [1 -1]; % AKA
AB1.s = [0 1];    % Euler

AB2.r = poly([0 1]);
AB2.s = [0 1.5 -0.5];

AB3.r = poly([0 0 1]);
AB3.s = [0 23/12 -4/3 5/12];

AB4.r = poly([0 0 0 1]);
AB4.s = [0 55/24 -59/24 37/24 -3/8];

AB5.r = poly([0 0 0 0 1]);
AB5.s = [0 1901/720 -1387/360 109/30 -637/360 251/720];

% Adams-Moulton

AM0.r=[1 -1]; % AKA
AM0.s=[1 0];  % backward Euler

AM1.r = [1 -1];    % AKA
AM1.s = [0.5 0.5]; % trapezoidal rule

AM2.r = poly([0 1]);
AM2.s = [5/12 2/3 -1/12];

AM3.r = poly([0 0 1]);
AM3.s = [3/8 19/24 -5/24 1/24];

AM4.r = poly([0 0 0 1]);
AM4.s = [251/720 646/720 -264/720 106/720 -19/720];

lmm = {BDF1,BDF2,BDF3, BDF4, BDF5, BDF6, AB1, AB2, AB3, AB4, AB5,...
              AM0, AM1, AM2, AM3, AM4 };
          
global os;

% ***One-step Methods***
% Defines the stability function R(z)

RK2E = '1+z+z^2/2';
RK2I = '(1+z/3)/(1-2*z/3+z^2/6)';
GL2 = '(1+z/2+z^2/12)/(1-z/2+z^2/12)';

os = {RK2E, RK2I, GL2};

% Last Modified by GUIDE v2.5 24-Aug-2012 11:24:38

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @stability_gui_OpeningFcn, ...
                   'gui_OutputFcn',  @stability_gui_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before stability_gui is made visible.
function stability_gui_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to stability_gui (see VARARGIN)

% Choose default command line output for stability_gui
handles.output = hObject;

axes(handles.axes3);
handles.rho_text=...
    text(0,0,'$\rho = $','Interpreter','latex','FontSize',16.0);

axes(handles.axes4);
handles.sigma_text=...
    text(0,0,'$\sigma = $','Interpreter','latex','FontSize',16.0);

axes(handles.axes5);
handles.R_text=...
    text(0,0,'$R(z) = $','Interpreter','latex','FontSize',16.0);
%set(handles.R_text,'Visible','off');

os=get(handles.select_one_step,'Value');
lmm=get(handles.select_lmm,'Value');

global check_axes
check_axes = 'custom';
% Update handles structure
guidata(hObject, handles);

% UIWAIT makes stability_gui wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = stability_gui_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;



function enter_rho_Callback(hObject, eventdata, handles)
% hObject    handle to enter_rho (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of enter_rho as text
%        str2double(get(hObject,'String')) returns contents of enter_rho as a double

% --- Executes during object creation, after setting all properties.
function enter_rho_CreateFcn(hObject, eventdata, handles)
% hObject    handle to enter_rho (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function enter_sigma_Callback(hObject, eventdata, handles)
% hObject    handle to enter_sigma (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of enter_sigma as text
%        str2double(get(hObject,'String')) returns contents of enter_sigma as a double

% --- Executes during object creation, after setting all properties.
function enter_sigma_CreateFcn(hObject, eventdata, handles)
% hObject    handle to enter_sigma (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on button press in select_lmm.
function select_lmm_Callback(hObject, eventdata, handles)
% hObject    handle to select_lmm (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of select_lmm

set(hObject,'Value',... toggle on this button
    get(hObject,'Max'));
set(handles.select_one_step,'Value',... set toggle off on other button
    get(handles.select_one_step,'Min')); 
    
set(handles.one_step_panel,'Visible','off');
set(handles.lmm_panel,'Visible','on');

% --- Executes on button press in select_one_step.
function select_one_step_Callback(hObject, eventdata, handles)
% hObject    handle to select_one_step (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of select_one_step

set(hObject,'Value',... toggle on this button
    get(hObject,'Max'));
set(handles.select_lmm,'Value',...
    get(handles.select_lmm,'Min')); % set toggle off on other button
set(handles.one_step_panel,'Visible','on');
set(handles.lmm_panel,'Visible','off');

    

function enter_R_Callback(hObject, eventdata, handles)
% hObject    handle to enter_R (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of enter_R as text
%        str2double(get(hObject,'String')) returns contents of enter_R as a double

% --- Executes during object creation, after setting all properties.
function enter_R_CreateFcn(hObject, eventdata, handles)
% hObject    handle to enter_R (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on button press in lmm_button.
function lmm_button_Callback(hObject, eventdata, handles)
% hObject    handle to lmm_button (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

try
    sigma=str2num(get(handles.enter_sigma,'String'));
    rho=str2num(get(handles.enter_rho,'String'));
    assert( all(size(rho)==size(sigma)) && ...
        and(size(rho,1)==1,size(sigma,1)) && ...
        and(isreal(rho),isreal(sigma)) && ...
        ~or(any(isnan(rho)),any(isnan(sigma))) &&...
        ~or(any(isinf(rho)),any(isinf(sigma))),...
        'Rho and sigma must be real row vectors of equal lengths!');

    % get x plotting range
    xlower=str2double(get(handles.xlower,'String'));
    xupper=str2double(get(handles.xupper,'String'));
    assert( ~or(isnan(xlower),isnan(xupper)) &&  xlower < xupper,...
    'Incorrect x plotting range!');

    %get y plotting range
    ylower=str2double(get(handles.ylower,'String'));
    yupper=str2double(get(handles.yupper,'String'));
    
    assert( ~or(isnan(ylower),isnan(yupper)) &&  ylower < yupper,...
    'Incorrect y plotting range!');

    % get resolution
    res=str2double(get(handles.res,'String'));
    assert( ~isnan(res) && isreal(res) && and(res <= 2000,res>=10) && ...
        floor(res) == res,...
    'Resolution must be an integer no less than 10 and no more than 2000');

    set(hObject,'String','Wait...');
    set(hObject,'Enable','off');
    drawnow;
    plot_lmm(handles.axes1,rho,sigma,[xlower xupper ylower yupper],res);
    set(hObject,'String','Plot');
    set(hObject,'Enable','on');
catch err
    cla(handles.axes1,'reset');
    h=msgbox(err.message,'Error','error','modal');
    waitfor(h);
end


% --- Executes on button press in one_step_button.
function one_step_button_Callback(hObject, eventdata, handles)
% hObject    handle to one_step_button (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

try
    % get x plotting range
    xlower=str2double(get(handles.xlower,'String'));
    xupper=str2double(get(handles.xupper,'String'));
    assert( ~or(isnan(xlower),isnan(xupper)) &&  xlower < xupper,...
    'Incorrect x plotting range!');

    %get y plotting range
    ylower=str2double(get(handles.ylower,'String'));
    yupper=str2double(get(handles.yupper,'String'));
    
    assert( ~or(isnan(ylower),isnan(yupper)) &&  ylower < yupper,...
    'Incorrect y plotting range!');

    % get resolution
    res=str2double(get(handles.res,'String'));
    assert( ~isnan(res) && isreal(res) && and(res <= 2000,res>=10) && ...
        floor(res) == res,...
    'Resolution must be an integer no less than 10 and no more than 2000');
    

    % get function R
    temp=inline(get(handles.enter_R,'String'),'z');
    temp=@(z)(temp(z));
    R=@(z)(arrayfun(temp,z));
    R=@(z)(R(z));

    R = inline(get(handles.enter_R,'string'),'z');
    z = sym('z');
    feval(R,z); % will trigger error if incorrect
    R=@(z)(R(z));
    R=@(z)(arrayfun(R,z));

    set(hObject,'String','Wait...');
    set(hObject,'Enable','off');
    drawnow;
    plot1step(handles.axes1,R,[xlower xupper ylower yupper],res);
    set(hObject,'String','Plot');
    set(hObject,'Enable','on');
catch err
    cla(handles.axes1,'reset');
    h=msgbox(err.message,'Error','error','modal');
    waitfor(h);
end


% --- Executes on selection change in popupmenu1.
function popupmenu1_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns popupmenu1 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu1
global lmm;

index = get(hObject,'Value');

if index ~= 1
    choice=lmm{index-1};
    
    % Reduce extra spaces to just two between values
    
    str=rats(choice.r);
    str=regexprep(str,'\s+','  ');
    set(handles.enter_rho,'String',['[' str ']']);
    
    str=rats(choice.s);
    str=regexprep(str,'\s+',' ');
    set(handles.enter_sigma,'String',['[' str ']']);
end

clear lmm;


% --- Executes during object creation, after setting all properties.
function popupmenu1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function xupper_Callback(hObject, eventdata, handles)
% hObject    handle to xupper (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of xupper as text
%        str2double(get(hObject,'String')) returns contents of xupper as a double



% --- Executes during object creation, after setting all properties.
function xupper_CreateFcn(hObject, eventdata, handles)
% hObject    handle to xupper (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function ylower_Callback(hObject, eventdata, handles)
% hObject    handle to ylower (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ylower as text
%        str2double(get(hObject,'String')) returns contents of ylower as a double



% --- Executes during object creation, after setting all properties.
function ylower_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ylower (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function yupper_Callback(hObject, eventdata, handles)
% hObject    handle to yupper (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of yupper as text
%        str2double(get(hObject,'String')) returns contents of yupper as a double



% --- Executes during object creation, after setting all properties.
function yupper_CreateFcn(hObject, eventdata, handles)
% hObject    handle to yupper (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function res_Callback(hObject, eventdata, handles)
% hObject    handle to res (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of res as text
%        str2double(get(hObject,'String')) returns contents of res as a double


% --- Executes during object creation, after setting all properties.
function res_CreateFcn(hObject, eventdata, handles)
% hObject    handle to res (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function xlower_Callback(hObject, eventdata, handles)
% hObject    handle to xlower (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of xlower as text
%        str2double(get(hObject,'String')) returns contents of xlower as a double



% --- Executes during object creation, after setting all properties.
function xlower_CreateFcn(hObject, eventdata, handles)
% hObject    handle to xlower (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on selection change in popupmenu2.
function popupmenu2_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns popupmenu2 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu2

global os;

index=get(hObject,'Value');

if(index ~=1)
    set(handles.enter_R,'String',os{index-1});
end

clear os;
    


% --- Executes during object creation, after setting all properties.
function popupmenu2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on button press in checkbox1.
function checkbox1_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox1
global check_axes
val = get(handles.checkbox1,'Value');

switch val
    case 1
        check_axes = 'default';
        set(handles.xlower,'enable','off');
        set(handles.ylower,'enable','off');
        set(handles.xupper,'enable','off');
        set(handles.yupper,'enable','off');
        
    case 0
        check_axes = 'custom';
        set(handles.xlower,'enable','on');
        set(handles.ylower,'enable','on');
        set(handles.xupper,'enable','on');
        set(handles.yupper,'enable','on');
end

function plot_lmm( ax, rho, sigma, lims, res  )
%plot_lmm Plot the linear stability domain (and root condition) for an LMM
%   Input arguments
%       ax = axes handle where plotting should occur
%       rho, sigma = define the coefficients of the scheme
%           These are respectively the LHS and RHS coefficients in
%           descending order
%       lims = plotting limits
%
% WARNING: May not work correctly if the boundary isn't a bounded curve
%
% The method used is based on page 162 (Section 7.6.1) of
%
% Finite Difference Methods for Ordinary and 
% Partial Differential Equations:
% Steady-State and Time-Dependent Problems
% by Randall J. LeVeque [SIAM, Philadelphia, 2007]

global check_axes

shade=[0 0 1]; % shade colour blue
empty=[1 1 1]; % empty colour white
boundary=[0 0 0]; % boundary colour black

clear i;
cla(ax);

% generate points on the unit circle by parametrisation
t = linspace(0,2*pi,res);
eit = exp(1i.*t);

% calculate points on the boundary of the stability domain
z = polyval(rho,eit)./polyval(sigma,eit);

x=real(z);
y=imag(z);

% plotting stuff
set(ax,'Color',empty); % set background colour; prevent weirdness
axes(ax);

% get all points of intersection of the curve with itself
[~,~,index]=intersections(x,y,1);
[~,~,back_index]=intersections(fliplr(x),fliplr(y),1);   
% second call is a workaround:
%   Only the index of the first intersection is returned, so we reverse the
%   arrays, effectively searching from the end to the beginning
back_index=length(z)-back_index;
index=[index;back_index];
index=sort(index); % put in increasing order

empty_region=true; % do we need to shade the background?


hold on;

if(~isempty(index)) % if there are any points of intersection
    
    for si=1:length(index) % do for each intersection in index
        
        next=mod(si,length(index))+1; % si+1, but with wraparound
        start=ceil(index(si));
        stop=floor(index(next));
        if(start<=stop) % range as usual
            vals=start:stop;
        else % range jumps over the end of the array back to the beginning
            vals=[1:stop start:length(z)];
        end
        %do plotting
        zz=mean(z(vals));
        p=rho-zz*sigma;
        
        if root_condition(p) % root condition
            % plot in blue (blue border); part of stability domain
            area(x(vals),y(vals),'FaceColor',shade,'EdgeColor','none');
            plot(x(vals),y(vals),'Color',boundary);
            empty_region=false;
        else
            % plot in white (blue border); not part of stability domain
            area(x(vals),y(vals),'FaceColor',empty,'EdgeColor','none');
            plot(x(vals),y(vals),'Color',boundary);
        end 
    end
else
    if( polyarea(x,y) < 10e-6 ) % if collinear points
        xmin=min(x);
        xmax=max(x);
        ymin=min(y);
        ymax=max(y);

        x0=(xmin+xmax)/2;
        y0=(ymin+ymax)/2;

        dx=xmax-xmin;
        dy=ymax-ymin;

        % normalise
        l=norm([dx dy]);

        dx=dx/l;
        dy=dy/l;

        % produce an orthogonal vector
        dxp=-dy;
        dyp=dx;

        % produce a trial point
        xt=x0+dxp;
        yt=y0+dyp;

        zt=xt+1i*yt;

        % check if trial point is in the stability domain
        p=rho-zt*sigma;

        if ~root_condition(p) % not root condition
        % go in opposite direction
            dxp=-dxp;
            dyp=-dyp;
        end

        % clip to fit within display window
        x1=min(lims(2),max(lims(1),xmax+l*dxp));
        x2=min(lims(2),max(lims(1),xmin+l*dxp));

        y1=min(lims(4),max(lims(3),ymax+l*dyp));
        y2=min(lims(4),max(lims(3),ymin+l*dyp));

        xmin=min(lims(2),max(lims(1),xmin));
        xmax=min(lims(2),max(lims(1),xmax));

        ymin=min(lims(4),max(lims(3),ymin));
        ymax=min(lims(4),max(lims(3),ymax));

        x=[xmin xmax x1 x2 xmin];
        y=[ymin ymax y1 y2 ymin];
        empty_region=false;
        fill(x,y,shade,'EdgeColor',boundary);
    else % have a region of reasonable area
        %do plotting
        zz=mean(z);
        %scatter(real(zz),imag(zz),'*g');
        p=rho-zz*sigma;

        if root_condition(p) % root condition
            % plot in blue (blue border); part of stability domain
            area(x,y,'FaceColor',shade,'EdgeColor','none');
            plot(x,y,'Color',boundary);
            empty_region=false;
        else
            % plot in white (blue border); not part of stability domain
            area(x,y,'FaceColor',empty,'EdgeColor','none');
            plot(x,y,'Color',boundary);
        end 
    end
end

if(empty_region) % is the stability domain everything outside the boundary?
   set(ax,'Color',shade);
end

refresh;

% Plot elements for visual inspection of root condition
plot(eit,'--m','LineWidth',2.0); % unit circle
r=roots(rho); % roots of rho
% following ensures correct plotting behaviour for all roots
if( r == real(r) ) % all real roots
    plot(r,zeros(size(r)),'or');
else % some imaginary roots
    plot(r,'or');
end


% axis settings 
switch check_axes
    case 'default' %default axes
        axis image; % equal x and y scales
        L = ylim;
        line([0 0],[L(1) L(2)],'Color','k','LineStyle','--','LineWidth',2.0);        
    case 'custom' %custom axes
        axis(lims)
        line([0 0],[lims(3) lims(4)],'Color','k','LineStyle','--','LineWidth',2.0);
end
hold off;


function cond = root_condition(p)
    rootsp=roots(p);
    % count of roots...
    in=0; % ...inside unit circle
    on=0; % ...on unit circle
    out=0; % ...outside unit circle

    for j=1:length(rootsp)
        if abs(rootsp(j)) > 1
            out=out+1;
        elseif abs(rootsp(j)) ==1
            on=on+1;
        else
            in=in+1;
        end
    end
    
    cond=(out == 0) && ~(on > 1);

    
function plot1step( ax, R, range, res )
%plot1step Plot the linear stability domain of certain 1 step methods
%   THEORY:
%   When applied to the test problem u'=bu, the typical 1 step method with
%   step size h will typically lead to the recurrence U(N+1)=R(z)U(N) for
%   the approximate solution U, where z=bh. 
%   
%   The method will be stable for exactly those z for which |R(z)| <= 1.
%
%   MATLAB's contourf function is used to plot filled contours. Some
%   adjustments are required to ensure that the correct region is shaded.
%   
%   INPUT:
%       ax ... handle to axes in which to do plotting
%       R ... a function handle, representing R(z) as described above
%       range ... represents the range of z to plot
%                 range = [xmin xmax ymin ymax];
%       resolution ... the number of points to use to partition the x range
%                      and y range.
%   OUTPUT:
%       h ... handle to contourgroup object produced
%
% The method used is based on page 163 (Section 7.6.2) of
%
% Finite Difference Methods for Ordinary and 
% Partial Differential Equations: Steady-State and Time-Dependent Problems
% by Randall J. LeVeque [SIAM, Philadelphia, 2007]

global check_axes

shade=[0 0 1]; % dark blue
empty=[1 1 1]; % white
boundary=[0 0 0];

[xx,yy]=meshgrid(linspace(range(1),range(2),res),...
    linspace(range(3),range(4),res));


% Generate values of |R(z)| for creating contours
zz=R(xx+1i*yy);
modzz=abs(zz);


[C,h]=contour(ax,xx,yy,modzz,[1 1]);
set(h,'LineColor',boundary); % make contour lines have 'boundary' colour


% This next bit extracts the x and y data for the contour from the contour
% matrix C.
%
% For more information on the contour matrix C: 
% 'doc contourc' and then click on the 'How To' link at the bottom.

index=1;
cc=1; % counts the number of contour 'records' in C

if(isempty(C))
    return;
end

% each record has a header of the contour level (1st row) and number of
% (x=1st row,y=2nd row) points which follow (2nd row).
%
% the latter value is the one we read into 'count'

count=C(2,index);

% read x and y from the record
currx=C(1,index+1:index+count);
curry=C(2,index+1:index+count);
% if there happen to be NaNs in currx or curry, throw them away and look
% for a new record (but not if we would go pas the end of the array
while (any(isnan(currx)) || any(isnan(curry)) ) && index+count+1<length(C)
    index=index+count+1; % move to next record
    cc=cc+1; % increase count of contour records
    count=C(2,index); % read new header
    currx=C(1,index+1:index+count); % get x from record
    curry=C(2,index+1:index+count); % get y from record
end

% if there are any NaNs left something has gone wrong and we shouldn't plot
if (any(isnan(currx)) || any(isnan(curry)) )
    return;
end

% the above code produces the x and y we want:
x=currx; 
y=curry;

% to shade, we have to test one point to see if it's inside or outside the 
% stability region. We choose the mean, since it will definitely lie
% on the inside of the contour if the contour is convex.
xbar=mean(x);
ybar=mean(y);

% work out the side on which each endpoint lies. Pattern:
%
%                     2
%       max y ---> ________
%                 |        |
%               1 |        | 3
%       min y --->|________|
%                  \   4    \
%                   \        \
%                    min x    max x

if(x(1) ~= x(end) || y(1) ~= y(end) ) % close off a non-closed contour
        perm=[1 3 4 2]; % permuting indices to get the numbering above
        inv_perm=[1 4 2 3]; % inverse
    
        tol=1e-3;

        sides_1=abs([x(1) x(1) y(1) y(1)]-range)<tol;
        side_1=find(sides_1);
        side_1=side_1(1:1); %(lying on multiple sides irrelevant)
        side_1=perm(side_1);

        sides_end=abs([x(end) x(end) y(end) y(end)]-range)<tol;
        side_end=find(sides_end);
        side_end=side_end(1:1); % take first one (lying on multiple sides irrelevant)
        side_end=perm(side_end);
            
        end0=length(x);
        switch( abs(side_1-side_end) )
            case 2 % opposite sides
                side_between=floor((side_1+side_end)/2);
                ind=inv_perm(side_between);
                
                if( ind <= 2 ) % determines x coordinate
                    x(end0+1)=range(ind);
                    y(end0+1)=y(end0);
                    x(end0+2)=range(ind);
                    y(end0+2)=y(1);
                    x(end0+3)=x(1);
                    y(end0+3)=y(1);
                else % determines y coordinate
                    x(end0+1)=x(end0);
                    y(end0+1)=range(ind);
                    x(end0+2)=x(1);
                    y(end0+2)=range(ind);
                    x(end0+3)=x(1);
                    y(end0+3)=y(1);
                end
                
            case 0 % same side
                x(end0+1)=x(1);
                y(end0+1)=y(1);
            otherwise % adjacent sides
                if side_end == 1 || side_end == 3
                    x(end0+1)=x(end0); % need to change y first (go along x)
                    y(end0+1)=y(1);
                    x(end0+2)=x(1);
                    y(end0+2)=y(1);
                else
                    x(end0+1)=x(1); % need to change x first (go along x)
                    y(end0+1)=y(end0);
                    x(end0+2)=x(1);
                    y(end0+2)=y(1);
                end
        end
end

xbar=mean(x);
ybar=mean(y);

if abs(R(xbar+1i*ybar))<=1
     fc=shade;
     bc=empty;
else
     fc=empty;
     bc=shade;
end

axes(ax);
hold on;
fill(x,y,fc,'EdgeColor',boundary);
hold off;

set(ax,'Color',bc);
drawnow;

xrange=max(abs(min(x)),abs(max(x)));
yrange=max(abs(min(y)),abs(max(y)));

% axis settings 
switch check_axes
    case 'default' %default axes
        axis image; % equal x and y scales 
        hold on;
        L = ylim;
        line([0 0],[L(1) L(2)],'Color','k','LineStyle','--','LineWidth',2.0);
        hold off;
    case 'custom' %custom axes
        axis( range )
        hold on;
        line([0 0],[range(3) range(4)],'Color','k','LineStyle','--','LineWidth',2.0);
        hold off;
end

% Copyright (c) 2008, Douglas M. Schwarz
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
% POSSIBILITY OF SUCH DAMAGE.
function [x0,y0,iout,jout] = intersections(x1,y1,x2,y2,robust)
%INTERSECTIONS Intersections of curves.
%   Computes the (x,y) locations where two curves intersect.  The curves
%   can be broken with NaNs or have vertical segments.
%
% Example:
%   [X0,Y0] = intersections(X1,Y1,X2,Y2,ROBUST);
%
% where X1 and Y1 are equal-length vectors of at least two points and
% represent curve 1.  Similarly, X2 and Y2 represent curve 2.
% X0 and Y0 are column vectors containing the points at which the two
% curves intersect.
%
% ROBUST (optional) set to 1 or true means to use a slight variation of the
% algorithm that might return duplicates of some intersection points, and
% then remove those duplicates.  The default is true, but since the
% algorithm is slightly slower you can set it to false if you know that
% your curves don't intersect at any segment boundaries.  Also, the robust
% version properly handles parallel and overlapping segments.
%
% The algorithm can return two additional vectors that indicate which
% segment pairs contain intersections and where they are:
%
%   [X0,Y0,I,J] = intersections(X1,Y1,X2,Y2,ROBUST);
%
% For each element of the vector I, I(k) = (segment number of (X1,Y1)) +
% (how far along this segment the intersection is).  For example, if I(k) =
% 45.25 then the intersection lies a quarter of the way between the line
% segment connecting (X1(45),Y1(45)) and (X1(46),Y1(46)).  Similarly for
% the vector J and the segments in (X2,Y2).
%
% You can also get intersections of a curve with itself.  Simply pass in
% only one curve, i.e.,
%
%   [X0,Y0] = intersections(X1,Y1,ROBUST);
%
% where, as before, ROBUST is optional.

% Version: 1.12, 27 January 2010
% Author:  Douglas M. Schwarz
% Email:   dmschwarz=ieee*org, dmschwarz=urgrad*rochester*edu
% Real_email = regexprep(Email,{'=','*'},{'@','.'})


% Theory of operation:
%
% Given two line segments, L1 and L2,
%
%   L1 endpoints:  (x1(1),y1(1)) and (x1(2),y1(2))
%   L2 endpoints:  (x2(1),y2(1)) and (x2(2),y2(2))
%
% we can write four equations with four unknowns and then solve them.  The
% four unknowns are t1, t2, x0 and y0, where (x0,y0) is the intersection of
% L1 and L2, t1 is the distance from the starting point of L1 to the
% intersection relative to the length of L1 and t2 is the distance from the
% starting point of L2 to the intersection relative to the length of L2.
%
% So, the four equations are
%
%    (x1(2) - x1(1))*t1 = x0 - x1(1)
%    (x2(2) - x2(1))*t2 = x0 - x2(1)
%    (y1(2) - y1(1))*t1 = y0 - y1(1)
%    (y2(2) - y2(1))*t2 = y0 - y2(1)
%
% Rearranging and writing in matrix form,
%
%  [x1(2)-x1(1)       0       -1   0;      [t1;      [-x1(1);
%        0       x2(2)-x2(1)  -1   0;   *   t2;   =   -x2(1);
%   y1(2)-y1(1)       0        0  -1;       x0;       -y1(1);
%        0       y2(2)-y2(1)   0  -1]       y0]       -y2(1)]
%
% Let's call that A*T = B.  We can solve for T with T = A\B.
%
% Once we have our solution we just have to look at t1 and t2 to determine
% whether L1 and L2 intersect.  If 0 <= t1 < 1 and 0 <= t2 < 1 then the two
% line segments cross and we can include (x0,y0) in the output.
%
% In principle, we have to perform this computation on every pair of line
% segments in the input data.  This can be quite a large number of pairs so
% we will reduce it by doing a simple preliminary check to eliminate line
% segment pairs that could not possibly cross.  The check is to look at the
% smallest enclosing rectangles (with sides parallel to the axes) for each
% line segment pair and see if they overlap.  If they do then we have to
% compute t1 and t2 (via the A\B computation) to see if the line segments
% cross, but if they don't then the line segments cannot cross.  In a
% typical application, this technique will eliminate most of the potential
% line segment pairs.


% Input checks.
error(nargchk(2,5,nargin))

% Adjustments when fewer than five arguments are supplied.
switch nargin
	case 2
		robust = true;
		x2 = x1;
		y2 = y1;
		self_intersect = true;
	case 3
		robust = x2;
		x2 = x1;
		y2 = y1;
		self_intersect = true;
	case 4
		robust = true;
		self_intersect = false;
	case 5
		self_intersect = false;
end

% x1 and y1 must be vectors with same number of points (at least 2).
if sum(size(x1) > 1) ~= 1 || sum(size(y1) > 1) ~= 1 || ...
		length(x1) ~= length(y1)
	error('X1 and Y1 must be equal-length vectors of at least 2 points.')
end
% x2 and y2 must be vectors with same number of points (at least 2).
if sum(size(x2) > 1) ~= 1 || sum(size(y2) > 1) ~= 1 || ...
		length(x2) ~= length(y2)
	error('X2 and Y2 must be equal-length vectors of at least 2 points.')
end


% Force all inputs to be column vectors.
x1 = x1(:);
y1 = y1(:);
x2 = x2(:);
y2 = y2(:);

% Compute number of line segments in each curve and some differences we'll
% need later.
n1 = length(x1) - 1;
n2 = length(x2) - 1;
xy1 = [x1 y1];
xy2 = [x2 y2];
dxy1 = diff(xy1);
dxy2 = diff(xy2);

% Determine the combinations of i and j where the rectangle enclosing the
% i'th line segment of curve 1 overlaps with the rectangle enclosing the
% j'th line segment of curve 2.
[i,j] = find(repmat(min(x1(1:end-1),x1(2:end)),1,n2) <= ...
	repmat(max(x2(1:end-1),x2(2:end)).',n1,1) & ...
	repmat(max(x1(1:end-1),x1(2:end)),1,n2) >= ...
	repmat(min(x2(1:end-1),x2(2:end)).',n1,1) & ...
	repmat(min(y1(1:end-1),y1(2:end)),1,n2) <= ...
	repmat(max(y2(1:end-1),y2(2:end)).',n1,1) & ...
	repmat(max(y1(1:end-1),y1(2:end)),1,n2) >= ...
	repmat(min(y2(1:end-1),y2(2:end)).',n1,1));

% Force i and j to be column vectors, even when their length is zero, i.e.,
% we want them to be 0-by-1 instead of 0-by-0.
i = reshape(i,[],1);
j = reshape(j,[],1);

% Find segments pairs which have at least one vertex = NaN and remove them.
% This line is a fast way of finding such segment pairs.  We take
% advantage of the fact that NaNs propagate through calculations, in
% particular subtraction (in the calculation of dxy1 and dxy2, which we
% need anyway) and addition.
% At the same time we can remove redundant combinations of i and j in the
% case of finding intersections of a line with itself.
if self_intersect
	remove = isnan(sum(dxy1(i,:) + dxy2(j,:),2)) | j <= i + 1;
else
	remove = isnan(sum(dxy1(i,:) + dxy2(j,:),2));
end
i(remove) = [];
j(remove) = [];

% Initialize matrices.  We'll put the T's and B's in matrices and use them
% one column at a time.  AA is a 3-D extension of A where we'll use one
% plane at a time.
n = length(i);
T = zeros(4,n);
AA = zeros(4,4,n);
AA([1 2],3,:) = -1;
AA([3 4],4,:) = -1;
AA([1 3],1,:) = dxy1(i,:).';
AA([2 4],2,:) = dxy2(j,:).';
B = -[x1(i) x2(j) y1(i) y2(j)].';

% Loop through possibilities.  Trap singularity warning and then use
% lastwarn to see if that plane of AA is near singular.  Process any such
% segment pairs to determine if they are colinear (overlap) or merely
% parallel.  That test consists of checking to see if one of the endpoints
% of the curve 2 segment lies on the curve 1 segment.  This is done by
% checking the cross product
%
%   (x1(2),y1(2)) - (x1(1),y1(1)) x (x2(2),y2(2)) - (x1(1),y1(1)).
%
% If this is close to zero then the segments overlap.

% If the robust option is false then we assume no two segment pairs are
% parallel and just go ahead and do the computation.  If A is ever singular
% a warning will appear.  This is faster and obviously you should use it
% only when you know you will never have overlapping or parallel segment
% pairs.

if robust
	overlap = false(n,1);
	warning_state = warning('off','MATLAB:singularMatrix');
	% Use try-catch to guarantee original warning state is restored.
	try
		lastwarn('')
		for k = 1:n
			T(:,k) = AA(:,:,k)\B(:,k);
			[unused,last_warn] = lastwarn;
			lastwarn('')
			if strcmp(last_warn,'MATLAB:singularMatrix')
				% Force in_range(k) to be false.
				T(1,k) = NaN;
				% Determine if these segments overlap or are just parallel.
				overlap(k) = rcond([dxy1(i(k),:);xy2(j(k),:) - xy1(i(k),:)]) < eps;
			end
		end
		warning(warning_state)
	catch err
		warning(warning_state)
		rethrow(err)
        return
	end
	% Find where t1 and t2 are between 0 and 1 and return the corresponding
	% x0 and y0 values.
	in_range = (T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) <= 1 & T(2,:) <= 1).';
	% For overlapping segment pairs the algorithm will return an
	% intersection point that is at the center of the overlapping region.
	if any(overlap)
		ia = i(overlap);
		ja = j(overlap);
		% set x0 and y0 to middle of overlapping region.
		T(3,overlap) = (max(min(x1(ia),x1(ia+1)),min(x2(ja),x2(ja+1))) + ...
			min(max(x1(ia),x1(ia+1)),max(x2(ja),x2(ja+1)))).'/2;
		T(4,overlap) = (max(min(y1(ia),y1(ia+1)),min(y2(ja),y2(ja+1))) + ...
			min(max(y1(ia),y1(ia+1)),max(y2(ja),y2(ja+1)))).'/2;
		selected = in_range | overlap;
	else
		selected = in_range;
	end
	xy0 = T(3:4,selected).';
	
	% Remove duplicate intersection points.
	[xy0,index] = unique(xy0,'rows');
	x0 = xy0(:,1);
	y0 = xy0(:,2);
	
	% Compute how far along each line segment the intersections are.
	if nargout > 2
		sel_index = find(selected);
		sel = sel_index(index);
		iout = i(sel) + T(1,sel).';
		jout = j(sel) + T(2,sel).';
	end
else % non-robust option
	for k = 1:n
		[L,U] = lu(AA(:,:,k));
		T(:,k) = U\(L\B(:,k));
	end
	
	% Find where t1 and t2 are between 0 and 1 and return the corresponding
	% x0 and y0 values.
    pos=(T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) < 1 & T(2,:) < 1).';
    neg=(T(1,:) <= 0 & T(2,:) <= 0 & T(1,:) > -1 & T(2,:) > -1).';
    
	in_range = (T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) < 1 & T(2,:) < 1).';
	x0 = T(3,in_range).';
	y0 = T(4,in_range).';
	
	% Compute how far along each line segment the intersections are.
	if nargout > 2
		iout = i(in_range) + T(1,in_range).';
		jout = j(in_range) + T(2,in_range).';
	end
end
