function F = force_beam(N, xstart, xend, beamlen, coeff, der);

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

h = (xend - xstart)./N; 
g = gauss_points(N, h, xstart);
f1 = gauss_weights(N, h)*coeff;

[bvn, bvn2] = bevaluate_beam(N, 4*N, g, beamlen./N);
if (der == 2)
    bvn = bvn2;
end

F = zeros(N+1, 1);
for n = 1:N+1
    F(n) = sum(f1.*bvn(:,n));
end
F = F;