How to improve uncertainties from fitting algorithm

I'm trying to fit some data (435 data points) to this equation:

$$\frac{1}{q^2}\sum_{k,l}^3\rho_k\rho_l\exp(-q^2(\sigma_k^2+\sigma_l^2)/2)\cos[q(d_k-d_l)] +bkg$$

where $q$ is the indipendent variable and everything else is a free parameter (with soft costraints in a range) for a total of 9 parameters. I perform the fitting with Matlab, using these options:

opts = fitoptions( 'Method', 'NonlinearLeastSquares' );

opts.Algorithm = 'Trust-Region';

opts.Robust = 'Bisquare';

opts.MaxIter = 1e8;

opts.MaxFunEvals = 1e8;

opts.TolFun = 1e-12;

opts.TolX = 1e-12;

opts.StartPoint = [0,0.1,-0.1,0.1,3,8,3,-18,18];

opts.Lower = [0,0,-0.1,0,0,0,0,-20,10];

opts.Upper = [0.003,0.1,0,0.1,8,15,8,-10,20];

what I get is a very good fit (visually, but also $R^2$ is very close to 1), but the 95% confidence bounds are extremely large, providing a relative error of $\sim 100\%$. Is there anything I can do to get better uncertainties? I suspect that the issue is that there are too many free parameters, but I don't know how to solve this problem.