carryOutNumericalIntegration = 1;    % set to 0 if not
meshpts = (0:1/51:1)';
%meshpts = [linspace(0,0.35,7)(1:6) linspace(0.35,0.65,200) linspace(0.65,1,7)(2:7)]';


function eff = forcingFn(x)
  eff = 100 * sin(sign(x - 0.5) .* exp(1 ./ (4*abs(x - 0.5).^1.05 + 0.3))) ...
          .* exp(1 ./ (4*abs(x - 0.5).^1.2 + 0.2) - 100*(x - 0.5).^2);
end

function M = stiffnessMat(meshPts)
  h_is = diff(meshPts);
  emm = diag( 1./h_is(1:end-1) + 1./h_is(2:end) );
  emm -= diag( 1./h_is(2:end-1), -1 );
  emm -= diag( 1./h_is(2:end-1), 1 );
  M = sparse(emm);
end

function bee = loadVecNumIntegrated(f, mPts)
  % load vector constructed one element at a time via numerical integration
  h_is = diff(mPts);
  for j=1:length(mPts) - 2
    bee(j) = quad(@(x) f(x)*(x - mPts(j)), mPts(j), mPts(j+1))/h_is(j) + ...
               quad(@(x) f(x)*(mPts(j+2) - x), mPts(j+1), mPts(j+2))/h_is(j+1);
  end
  bee = bee';
end

function bee = loadVecTrapRule(f, meshPts)
  % load vector constructed all-at-once via Trapezoid rule approximation
  h_is = diff(meshPts);
  bee = f(meshPts(2:end-1)) .* (h_is(1:end-1) + h_is(2:end)) / 2;
end

if carryOutNumericalIntegration
  b = loadVecNumIntegrated(@forcingFn, meshpts);
else
  b = loadVecTrapRule(@forcingFn, meshpts);
end
plot(meshpts, [0; stiffnessMat(meshpts) \ b; 0])