how modify the levenberg Marquardt parameter(lambda0)

hello,

my levenberg Marquardt code to identify 2 parameter is showing below:

normally levenberg Marquardt method acts like the Steepest Descent method when the parameters are far from their optimal value, and acts like the Gauss Newton method when the parameters are close to the optimal value. but in my example i have not seen that lavenberg marquardt's method varies between Steepest Descent method and Gauss Newton method when i change the initial value of two variables and the initial value of lambda0.

maybe there is a parameter to rule, can you check the code with me please

function LevenbergMarquardtMethod %p est un vecteur de dimension 2x1 contenant les parametres
              clear all; clc;
              p(1,1) = 0.0004*200;
              p(2,1) = 60+200;
              %F=Rexp-p(1)*(1-exp(-t/p(2))); -> forme de l'equation a minimiser 
              %dF1 = -(1-exp(-t/p(2))); -> forme de la premiere derivee en p(1)
              %dF2 = (p(1)*t/(p(2)^2)).*(exp(-t/p(2))); -> forme de la seconde derivee en p(2)

                Rexact=0.0004;
                tauxexact=60;
                t = [ 0; 2; 4; 6; 8; 10; 12; 14; 16; 18; 20; 22; 24; 26; 28; 30; 32;34;36;38;40;42;44;46;48;50;52;54;56;58;60;62;64;66;68;70;72;74;76;78;80;82;84;86;88;90;92;94;96;98;100;102;104;106;108;110;112;114;116;118;120];
                Rexp=Rexact*(1-exp(-t/tauxexact));
                %initialisation de lambda 
                lambda0 =0.01;
                rmsea = 15;
                rmse = 10;
                inc = 0;
                %for inc = 1:100
                while abs(rmsea-rmse(end))>1e-5 && inc<500 && lambda0 < 1e7 
                    inc = inc+1;
                    %calcul du jacobien 
                    J(1:length(t),1) = -(1-exp(-t/p(2)));
                    J(1:length(t),2) = (p(1)*t/(p(2)^2)).*(exp(-t/p(2)));
                    %calcul des residus
                    F=Rexp-p(1)*(1-exp(-t/p(2)));
                    %calul de la RMSE lors de la premiere iteration 
                    if inc==1  
                        rmse = (sum(F.^2)/length(F))^0.5;
                    else
                        %mise a jour de lambda 
                        lambda = lambda0
                        %generation de la matrice unite 
                        tailleJJ = length(J.'*J);

                        %calul des delta en fonction de lambda 
                        Delta1 = -inv((J.'*J)+lambda*eye(tailleJJ))*(J.')*F;
                        lambda=lambda0/10;
                        Delta2 = -inv((J.'*J)+lambda*eye(tailleJJ))*(J.')*F;
                        p1 = p; 
                        p2 = p;
                        %mise a jour des parametres dans un nouveau vecteur 
                        p1 = p1+Delta1;
                        %calcul des nouveaux residus et de la nouvelle 
                        %RMSE pour lambda=lambda0
                        F=Rexp-p1(1)*(1-exp(-t/p1(2)));
                        rmse1 = (sum(F.^2)/length(F))^0.5; 
                        s=sum(F.^2)
                        %calcul des nouveaux residus et de la nouvelle %RMSE pour lambda=lambda0/10
                        p2 = p2+Delta2;
                F=Rexp-p2(1)*(1-exp(-t/p2(2))); 
                rmse2 = (sum(F.^2)/length(F))^0.5;
                %condition sur la mise a jour des parametres
                %voir texte pour plus de details 
                if rmse1>rmse && rmse2>rmse 

                lambda0 = lambda0*10;
                    elseif rmse2<rmse1

                        lambda0 = lambda0/10;
                        p = p2;
                        rmsea = rmse; 
                        rmse = rmse2;
                    elseif rmse2>=rmse1

                        p = p1; 
                        rmsea = rmse ;
                        rmse = rmse1;
                    end
                    end

                end

                %calcul de la RMSE
                rmse = (sum(F.^2)/length(F))^0.5;
                %calcul de l'ecart type (standard deviation) 
                stdv = (sum(F.^2)/(length(F)-length(p)))^0.5;
                %variance
                segma=((sum((F.^2)/rmse))/(length(F)-length(p)))^0.5;
                r=(((60-p(2))^2+(0.0004-p(1))^2)/2)^0.5;
                %evalutation des donnees ajustees 
                ta = 0:0.01:120;
                %Da = (p(1).*ta)./(p(2)+ta);
                Da = p(1)*(1-exp(-ta/p(2)));

                %affichage des resultats 
                R=p(1)
                taux=p(2)
                %fprintf('\np1 = %f\ntp2 = %f\n\n',p);
                erreur_R=(abs(Rexact-p(1))/Rexact)*100
                erreur_taux=(abs(tauxexact-p(2))/tauxexact)*100
                fprintf('rmse = %f\n',rmse); 
                fprintf('stdv = %f\n\n',stdv);
                fprintf('segma = %f\n\n',segma);
                fprintf('résidu = %f\n\n',r);
                %affichage du graphique figure
                h = plot(t, Rexp, ':k*', ta, Da, 'r'); 
                legend(h, {'Donnees experimentales' 'Meilleur ajustement'}, 'location', 'southeast') 
                title('Ajustement avec la methode de Levenberg-Marquardt'); 
                xlabel('t'); 
                ylabel('Donnees a ajuster');