function [m,npath,eigval] = solve(m,varargin)
% solve  Calculate first-order accurate solution of the model.
%
% Syntax
% =======
%
%     M = solve(M,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Paramterised model object. Non-linear models must also
% have a steady state values assigned.
%
% Output arguments
% =================
%
% * `M` [ model ] - Model with newly computed solution.
%
% Options
% ========
%
% * `'expand='` [ numeric | *`0`* ] - Number of periods ahead up to which the
% model solution will be expanded.
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar in the command
% window.
%
% * `'refresh='` [ *`true`* | `false` ] - Refresh dynamic links before computing
% the solution.
%
% * `'select='` [ *`true`* | `false` ] - Automatically detect which equations
% need to be re-differentiated based on parameter changes from the last
% time the model was solved.
%
% * `'warning='` [ *`true`* | `false` ] - Display warnings produced by this
% function.
%
% Description
% ============
%
% The IRIS solver uses an ordered QZ (or generalised Schur) decomposition
% to integrate out future expectations. The QZ may (very rarely) fail for
% numerical reasons. IRIS  includes two patches to handle the some of the
% QZ failures: a SEVN2 patch (Sum-of-EigenValues-Near-Two), and an E2C2S
% patch (Eigenvalues-Too-Close-To-Swap).
%
% * The SEVN2 patch: The model contains two or more unit roots, and the QZ
% algorithm interprets some of them incorrectly as pairs of eigenvalues
% that sum up accurately to 2, but with one of them significantly below 1
% and the other significantly above 1. IRIS replaces the entries on the
% diagonal of one of the QZ factor matrices with numbers that evaluate to
% two unit roots.
%
% * The E2C2S patch: The re-ordering of thq QZ matrices fails with a
% warning `'Reordering failed because some eigenvalues are too close to
% swap.'` IRIS attempts to re-order the equations until QZ works. The
% number of attempts is limited to `N-1` at most where `N` is the total
% number of equations.
%
% Example
% ========
%

% -IRIS Toolbox. 2008/10/20.
% -Copyright 2007-2012 Jaromir Benes.

opt = passvalopt('model.solve',varargin{1:end});

% Backward compatibility.
if opt.forward > 0 && opt.expand == 0
    opt.expand = opt.forward;
end

%**************************************************************************

% Refresh dynamic links.
if opt.refresh && ~isempty(m.Refresh)
    m = refresh(m);
end

if opt.warning
    % Warning if some log-lin variables have non-positive steady state.
    chk(m,Inf,'parameters','log');
end

% Find solutions for all parameterisations.
[m,npath,sing2] = mysolve(m,Inf,opt);

if (opt.warning || opt.error) && any(npath ~= 1) && iswarning('model')
    % NaN eigenvalues.
    if opt.error
        func = @utils.error;
    else
        func = @utils.warning;
    end
    printfunc = @(x) sprintf(' #%g',find(x));
    index = npath == -2;
    if any(index)
        func('model', ...
            'Solution not available. Linear dependency in some equations in%s.', ...
            printfunc(index));
    end
    % No stable solution.
    index = npath == 0;
    if any(index)
        func('model', ...
            'Solution not available. No stable solution in%s.', ...
            printfunc(index));
    end
    % Infinitely many stable solutions.
    index = isinf(npath);
    if any(index)
        func('model', ...
            'Solution not available. Multiple stable solutions in%s.', ...
            printfunc(index));
    end
    % Complex derivatives.
    index = imag(npath) ~= 0;
    if any(index)
        func('model', ...
            'Solution not available. Complex derivatives in%s.', ...
            printfunc(index));
    end
    % NaN derivatives.
    index = isnan(npath);
    if any(index)
        func('model', ...
            'Solution not available. NaN derivatives in%s.', ...
            printfunc(index));
    end
    % Singularity in state space.
    index = npath == -1;
    if any(index)
        if any(sing2(:))
            dosingmeasurement();
        else
            func('model', ...
                ['Solution not available. ', ...
                'Singularity in state-space matrices in%s.'], ...
                printfunc(index));
        end
    end
end

if nargout > 2
    eigval = m.eigval;
end

% Nested functions.

%**************************************************************************
    function dosingmeasurement()
        pos = find(any(sing2,2));
        pos = pos(:).';
        report = {};
        for ieq = pos
            report{end+1} = strfun.alt2str(sing2(ieq,:)); %#ok<AGROW>
            report{end+1} = m.eqtn{ieq}; %#ok<AGROW>
        end
        func('model', ...
            ['Solution not available. ', ...
            'Singularity or NaN in this measurement equation in%s: ''%s''.'], ...
            report{:});
    end
% dosingmeasurement().

end