function Q = qmatrix(t, k)
% QMATRIX   Matrix of coefficients of the trapezoidal rule.
%    Q = QMATRIX(N), where N is a scalar, returns an NxN matrix whose
%    elements are the coefficients of the trapezoidal rule corresponding
%    to equidistant abscissa.
% 
%    Q = QMATRIX(T), where T is an N-dimensional vector, returns an NxN
%    matrix whose elements are the coefficients of the trapezoidal rule
%    corresponding to the abscissa specified by T. The elements of T
%    do not need to be equidistant.
%    
%    The cumulative integral of a vector X containing values correspoding 
%    to abscissa stored in a vector T can be computed as follows:
%    I = QMATRIX(T) * X
% 
%    Q = QMATRIX(N, K) or Q = QMATRIX(T, K) generates analogous matrices
%    but representative of the trapezoidal rule applied K times.

if ~isscalar(k) || floor(k) ~= k
    error('k must be a scalar integer number.')
end

if ~exist('k', 'var'), k = 1; end   % default: single integration

if isscalar(t), t = (1:t)'; end     % number of points specified instead
                                    % of vector of abscissa
N = length(t);
    
if k == 0, Q = eye(N);              % trapezoidal rule applied 0 times
                                    % => identity matrix
elseif k >= 1
    
    % Coefficients of the trapezoidal rule
    Q = zeros(N);
    for n = 2:N
        Q(n:N, [n-1, n]) = Q(n:N, [n-1, n]) + (t(n)-t(n-1))/2;
    end
    
    % Trapezoidal rule applied k times
    Q = Q^k;
end