The Nelson-Siegel function is non-linear. However, if you know the value of tau, then it is possible to estimate the betas with linear regression. That is what the function lsbetas does. You specify tau, then lsbetas performas linear regression and returns two variables: the coefficients (betas) and the sum of the squared residuals (e' * e). The first line of the function nelsonfit calls an optimizing function that searches for the value of tau that minimizes the sum of the squared residuals. The smaller the latter the better the fit.
You could also use the function lsqnonlin to estimate the betas and the tau. As the English say, there's more than one way to skin a cat.