How does lsqcurvefit estimate error at each iteration of parameter estimation?
It both checks how far it has moved (i.e. am I not moving very much anymore? Then stop), and the norm of the gradient (minima are when the gradient is zero). I think there's a third tolerance, "function tolerance", which I don't quite know. It's all here in the docs (search for tolerance).
Why does lsqcurvefit not find the same solution each time?
Levenberg-Marquardt is a local nonlinear least squares optimizer, so it's dependent on the initial condition (hence the term local). The reason is because it uses a Newton method (which is clearly local: it falls to the closest minima) mixed with gradient decent (it always goes downhill, hence going to the local minima). The result it that it finds minima closest to where you start, so it has initial condition dependence.
What approaches should be taken to choose a solution?
You should take a look at global optimization algorithms. MATLAB's Global Optimization Toolbox is okay, but I find the litany of methods provided by Julia's JuMP and Optim.jl to be better suited for most problems. These give you access to a whole range of packages, and you can use the same code to call all of the different solvers. I have tended to have most success with some methods from NLopt, though IPOPT can do well too.