function [s0,A,rho,rhos]=max_eig_matrix(Afun,s,order,maxiter,tol)
% MAX_EIG_MATRIX
%
% A = max_eig_matrix(Afun,s,order)
%

if ~exist('maxiter','var') || isempty(maxiter), maxiter=100; end
if ~exist('tol','var') || isempty(tol), tol=1e-6; end

assert(order>2);

p = cell(size(s));
n = order;
if isscalar(n), n=n*ones(size(s)); end
for i=1:numel(s)
    p{i} = gaussian_quadrature(s(i).recur(n(i)),n(i));
    p{i} = [s(i).l; p{i}; s(i).r];
end

if numel(p)==1,
    % dereference the cell if there is only one element
    p=p{1};
else
    % build a tensor grid of points
    pgrid=1; 
    for i=1:numel(p)
        pgrid=[kron(pgrid,ones(length(p{i}),1)) kron(ones(size(pgrid,1),1),p{i})];
    end
    p=pgrid(:,2:end);
end

% we need to be slightly ineffeicient to make this parfor friendly
rhos = zeros(size(p,1),1);
parfor i=1:size(p,1)
    rhos(i) = maxeig(Afun(p(i,:)),maxiter,tol);
end

[rho, maxind] = max(rhos);
s0 = p(maxind,:);
A = Afun(s0);


