@PeterSimon, @certik, thanks for the demonstration of the custom weight calculation. This is one of the big advantages of Julia, that you can leverage custom types across the whole library ecosystem.
As you’ve both noticed, the function f(x) (or M(x)) grows linearly as x approaches infinity. It’s supposed to model the probability of cell division. Above the critical mass x_c the function slopes upward.

However the model cell distribution behaves like exp(-x) for large x values (large cells are very rare, because they divide), so the weight in the integral actually looks something like this (depending on the trial function and weighting functions used in the method of weighted residuals):

My current plan was to first compare the results of Gauss-Laguerre quadrature with some adaptive methods available in QUADPACK, and MATLAB, to determine the accuracy. If not sufficient then I know custom quadrature rules are an option.