lsqnonlin nonlinear optimization and looping the initial parameters?

I am having trouble looping with a least squares nonlinear optimization routine. I am trying to fit some experimental data to my analytic function that I have created. The routine itself works fine (hsquared is consecutively added at updated times as a running sum, and each time the updated value is stored into the next row of a 360x1 vector), but the problem lies in the initial parameters. Since I am not sure what the parameters should be close to, I randomized them, but often I would encounter a local minimum and so I cannot just change them arbitrarily by small increments as this would take too long. I want to put the initial parameters in a loop so that each element of the initial parameter vector continues changing by set values, systematically, until the condition that I have specified is not longer met: that is, DiffOut > 0.15e7 (0.15e7 is the error threshold, and DiffOut is the vector of differences between the experimental data and the analytic function evaluated at each of those times from the data). In other words, once DiffOut becomes less than that value, I want to store the set of initial parameters params0 that caused that vector of differences to be less than 0.15e7, and use those for the remaining calculations. Note that DiffOut is only defined in the function file (correct?) Can anyone check over my code to see why this isn't working the way I am describing it? One file is a function .m file and the other is a script .m file. I am not sure if my looping structure is correct, and where to place the loop for the initial parameters (inside the other loop with the times?). I tried using a while statement since I want it to continue indefinitely until the condition is no longer satisfied, but Matlab does not allow while-else conditions so I tried using a for loop before the while, and putting an else after the while statement is closed. If you want the data file also, I can send that over. Please help!!

Here is the script:

clc clear

% Import the data sets that you are trying to fit the function to. Text file must be in same directory as this script

load d1.txt;

t=d1(:,1); h=d1(:,2); F=d1(:,3);

% Initialize params for two Voigt elements in series with one dashpot

params0=[rand(1) rand(1) rand(1) rand(1) rand(1)]';

% Set an options file for LSQNONLIN to use the medium-scale algorithm

options = optimset('Algorithm','levenberg-marquardt','TolX',1e-9,'MaxFunEvals',1200);

% Calculate the new coefficients using LSQNONLIN

params=lsqnonlin(@fit_discrete,params0,[1e-16 1e-16 1e-16 1e-16 1e-16]',[],options,t,h);

DiffOut = fit_discrete(params,t,h);

            % while each row of the difference vector is greater than the
            % error tolerance, iterate the initial conditions.  Otherwise,
            % once the condition is no longer satisfied and DiffOut <
            % 0.5e7 as required for the optimization, use the params
            % that led to that condition and set those as params_new

for k=1:50:300

for l=1:50:300

     for m=1:50:300

        if DiffOut > 0.5e7

            while DiffOut > 0.5e7

                DiffOut = fit_discrete(params,t,h);

                params_new=[params0(1)*k params0(2)*l params0(3)*m];

            end

            else 

                params_new=params;

            end

        end

    end
end

hsquared_new=zeros(size(t)); % initialize hsquared_new

for i=1:359

   for j=1:359

       if t(j) >= t(i)

             % creep compliance model based on dashpot connected in series to two Voigt elements

             hsquared_new(j)=hsquared_new(j)+(pi/4)*((F(i+1)-F(i))*((t(j)-...
                 t(i))/params_new(1)+(1/params_new(3))*(1-exp(-(t(j)-t(i))*params_new(3)/params_new(2)))+...
                 (1/params_new(5))*(1-exp(-(t(j)-t(i))*params_new(5)/params_new(4)))));

       else

             hsquared_new(j)=hsquared_new(j);

       end

   end

end

    % Plot the original and experimental data
            plot(t,h.^2,'+r',t,hsquared_new,'b')
            xlabel('t (s)');
            ylabel('h_c^2 (nm^2)');
            axis([0 30.2 0 2.5e7])

eta0_minimized=params(1) eta1_minimized=params(2) E1_minimized=params(3) eta2_minimized=params(4) E2_minimized=params(5)

And here is the function .m file.

function DiffOut = fit_discrete(params,t,h)

% This function is called by LSQNONLIN. % params is a vector which contains the coefficients of the equations. t % and h are the data sets that were passed to lsqnonlin.

eta0=params(1); eta1=params(2); E1=params(3); eta2=params(4); E2=params(5);

load d1.txt;

t=d1(:,1); % units of s h=d1(:,2); % units of nm F=d1(:,3); % units of uN

% creep compliance model based on dashpot connected in series to two Voigt % elements

hsquared=zeros(size(t));

for i=1:359

   for j=1:359

       if t(j) >= t(i)

           hsquared(j)=hsquared(j)+(pi/4)*((F(i+1)-F(i))*((t(j)-t(i))/eta0+(1/E1)*(1-exp(-(t(j)-t(i))*E1/eta1))+(1/E2)*(1-exp(-(t(j)-t(i))*E2/eta2))));

       else

           hsquared(j)=hsquared(j);

       end

   end

end

 DiffOut=hsquared-h.^2

end

[Star Strider]

Maybe I'm missing something here, but if you're fitting hsquared to (h.^2), why not use lsqcurvefit?

[Walter Roberson]

Please format your code. http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

[Hugo Luis Manterola]

http://www.mathworks.com/help/gads/multistart-using-lsqcurvefit-or-lsqnonlin.html check this out. I guess it's what you need. I also think you should use lsqcurvefit.