function [x0,y0] = ifplot(fun,lims,varargin)
%__________________________________________________________________________
%IFPLOT     Plot function specified to interval equations
%--------------------------------------------------------------------------
% * IFPLOT(FUN,LIMS) plots the interval function FUN between the x-axis
%   limits by LIMS=[XMIN XMAX], and the range estimation of FUN over
%   interval [XMIN,XMAX] as the y-axis limits. Using LIMS = [XMIN XMAX YMIN
%   YMAX] also controls the y-axis limits but without range estimation. The
%   function FUN(x) must return an interval row vector for each element of
%   vector x. For example, if FUN returns [f1(x), f2(x), f3(x)] then for
%   input [x1;x2] the function should return the matrix
%           
%           f1(x1)  f2(x1)  f3(x1)
%           f1(x2)  f2(x2)  f3(x2)
%
%   In case of real function FUN, FPLOT matlab function is utilized.
% -------------------------------------------------------------------------
%   Ioannis A. Nikas October 20, 2006
%   Update........................... November 23, 2007
%__________________________________________________________________________
%
error(nargchk(2,nargin,nargin)) % Validation of number of input arguments
Xo = infsup(min(lims(1:2)) , max(lims(1:2)));    % Define Xo
Yo = feval(fun,mid(Xo));                         % F(Xo)
marker = '';tol = 2e-3;N = 300;
%
if nargin >= 3 & ~isempty(varargin{1})
    if isstr(varargin{1})
        marker = varargin{1};
    elseif varargin{1} < 1
        tol = varargin{1};
    else
        N = varargin{1};
    end
end
if nargin >= 4 & ~isempty(varargin{2})
    if isstr(varargin{2})
        marker = varargin{2};
    elseif varargin{2} < 1
        tol = varargin{2};
    else 
        N = varargin{2};
    end
end
if nargin >= 5 & ~isempty(varargin{3})
    if isstr(varargin{3})
        marker = varargin{3};
    elseif varargin{3} < 1
        tol = varargin{3};
    else
        N = varargin{3};
    end
end
%
if ~isa(Yo,'intval')        % Real FUNCTION
    if any([strfind(marker,'|'), strfind(marker,'::'),strfind(marker,'=')]), error(msg(1)); end;
    fprintf(msg(2)); fplot(fun,lims,varargin{:}); box off;
else
    if (isvarname(fun))
        fun = ezfcnchk(fun);
    else
        fun = fcnchk(fun);
    end
    %
    % compute the x duration and minimum and maximum x step
    xmin = Xo.inf; xmax = Xo.sup;
    maxstep = diam(Xo)/N; minstep = min(maxstep, diam(Xo)*tol);
    tryVal = minstep;
    % compute the first two points
    x = xmin; y = feval(fun,x,varargin{4:end});
    xx = x;
    x = xmin+minstep;
    y(2,:) = feval(fun,x,varargin{4:end});
    xx(2) = x;
    % define y-axis limits either
    if length(lims) == 4    % by user
        ymin = min(lims(3:4)); ymax = max(lims(3:4));
        ylims = 1;
    else                    % or by range estimation
        Yrg = feval(fun,Xo);
        ymin = min(Yrg.inf); ymax = max(Yrg.sup);
        ylims = 0;
    end
    %
    I = 2;
    while xx(I) < xmax
        I = I+1;
        tryVal = min(maxstep,min(2*tryVal, xmax-xx(I-1)));  % Update Step
        x = xx(I-1) + tryVal;                               % New x
        y(I,:) = feval(fun,x,varargin{4:end});              % New y
        xx(I) = x;
    end
    %    
    if nargout == 0
        box on
        iplot(xx',y,marker);
        set(gca,'XLim',[xmin xmax]);
        set(gca,'YLim',[ymin ymax]);
    else
        x0 = xx; y0 = y;
    end
end
hold on
step = diam(Xo)*0.025;
stepf = diam(fun(Xo))*0.0125;
mn = inf(fun(mid(Xo)));
mx = sup(fun(mid(Xo)));
for x = Xo.inf:step:Xo.sup-step
    X = infsup(x,x+step);
    f = fun(X);
    mn = min(mn,f.inf);
    mx = max(mx,f.sup);
    line([X.inf,X.sup],[f.inf, f.inf],'Color',[.6 .4 .3])
    line([X.inf,X.sup],[f.sup, f.sup],'Color',[.6 .4 .3])
    line([X.inf,X.inf],[f.inf, f.sup],'Color',[.6 .4 .3])
    line([X.sup,X.sup],[f.inf, f.sup],'Color',[.6 .4 .3])
    drawnow
end

step = 2*step;
stepf = 2*stepf;
line([Xo.inf-step, Xo.sup + step],[0,0],'Color',[.6 .4 .3])
line([0, 0],[mn - stepf, mx + stepf],'Color',[.6 .4 .3])
axis([Xo.inf-step Xo.sup + step mn - stepf mx + stepf])

%__________________________________________________________________________
%==========================================================================
%--------------------------------------------------------------------------

function out = msg(error_code)

error(nargchk(1,nargin,nargin))
%
y(1) = {'IFPLOT ERROR: ''|'', ''|:'' and ''|:'' line styles are not supported by FPLOT function\n'};
y(2) = {'IFPLOT WARNING: FPLOT matlab function is utilized.\n'};
%
out = char(y(error_code));