The lambertw function could be a problem with respect to fitting your equation system. It is not always possible to use fminsearch (even though it is ‘derivative-free’) and other such optimisation functions and achieve a reasonable result.
If you have the Global Optimization Toolbox, I have found the ga function (see this Comment to HELP ME PLEASE! How to find parameter estimation with fminsearch? (link)) to provide good results. I used a very large, random, initial population. I got a very good fit to the data. It took a while, and as I see it, the time was well spent.
Gradient-descent nonlinear parameter estimation routines are sensitive to the initial parameter estimates, whereas ga, given a large initial population and enough time, will converge on the best result possible. If you don’t have access to the ga function, simple unconstrained genetic algorithm codes are relatively easy to write. (If it is not already obvious, I really like genetic algorithm optimisations!)
I have not run your code, or experimented with using ga to optimise it. I leave this for you to experiment with first.
