function A=mean_matrix(Afun,parameters,order)
% MEAN_MATRIX Output the average of a parameterized matrix.
%
% A = mean_matrix(Afun, parameters, order)
% return A which approximates the result of integrating Afun over 
% the parameter space using an "order"-order quadrature rule.

[p,w]=gaussian_quadrature(parameters,order);
A=Afun(p(1,:))*w(1);
parfor i=2:size(p,1)
    A = A + Afun(p(i,:))*w(i);
end
