For calculating the Liu and Thorn U_0 + DU, you need to sum over the first N zeros of the firstorder Bessel function. Unfortunately, Matlab doesn't seem to come with a function to do this.
Rather than reinvent the wheel, we can just use the function used by GWINC. 
function x=besselzero(n,k,kind)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% besselzero.m
%
% Find first k positive zeros of the Bessel function J(n,x) or Y(n,x)
% using Halley's method.
%
% Written by: Greg von Winckel  01/25/05
... 79 more lines ...
