function [M,K,C] = matrix_construct_beam(N, xstart, xend, beamlen, rho, EI, cD, gamma);

% This file is a helper to beam_model_demo.  Please see there for
% documentation

h = (xend - xstart)/N;
g = gauss_points(N, h, xstart);
[bvn, bvn2] = bevaluate_beam(N, 4*N, g, beamlen./N); % Note bvn2 is the 2nd derivative of bvn
wt  = gauss_weights(N, h);

Spline0 = zeros(N+1, N+1);
Spline2 = zeros(N+1, N+1);
for n=1:N+1
    for q = 1:N+1
        Spline0(q,n) = sum(wt .* bvn(:, q) .* bvn(:, n));
        Spline2(q,n) = sum(wt .* bvn2(:, q) .* bvn2(:, n));
    end
end

M = rho*Spline0;
K = EI*Spline2;
C = cD*Spline2 + gamma*Spline0;