I suggest that you learn to use the debugger. I put a break point in the solve_nch function, and found that at the first function call the value of A1 was about 1e-17 and nch was a vector of 1. Therefore exp(A1*nch) = 1 and so the logarithm evaluated to -Inf.
Alan Weiss
MATLAB mathematical toolbox documentation
P.S. Perhaps I should point out that for small a, , so