Optimizing/Fitting 2 variables

Hi everyone!

So far I've optimized the following code for 1 parameter only : K_dep. It workds great.

function K_dep_fit=onevar
options = optimset('PlotFcns',@optimplotfval); 
K_dep_fit = fminsearch(@criterion, 0.0964, options); %initial value for K_dep=0.0964  
   function crit = criterion (K_dep)
     
       %%%%%constants%%%%%   
y_exp=[7.595131579, 8.238253968,8.236679198,8.582412281,9.562149123]; 
x_mlad1=1.2e-2;x_mlad2=1.7e-2;x_mlad3=2.2e-2;x_mlad4=2.7e-2;x_mlad5=3.2e-2;x_mlad=[x_mlad1, x_mlad2, x_mlad3, x_mlad4, x_mlad5]; 
V_exp=100; 
V=V_exp/60*10e-9; 
area = 289632*1e-12; OD =0.02; C0 = OD * 2.56e14  ;a2 = 0.52e-6; mu=8.9e-4;
g=9.81; deltap=51.793;Vs=2/9*deltap*g*a2*a2/mu; Vs2=6.003e-8; w=0.5e-2;
h=0.03e-2;l=5.3e-2;sigma=6*V/(w*h*h); sigma2=180.84;Kb=1.38064852e-23; T=293.15; 
D=Kb*T/(6*pi*a2*mu); D2=3.897e-13; 

k1=K_dep*D/a2; 
beta1=Vs./k1; 
gamma1=beta1-1; 
x_tilde=x_mlad.*Vs.^3./(sigma.*D.^2);
y_exp_dim=y_exp./area./80./C0./Vs;  

    y_calc=getflux(x_tilde,gamma1); 
    crit=sum(((y_exp_dim-y_calc)./y_exp_dim).^2);    

end

x_tilde5=0.00001:0.01:4;

y_calc_plot=getflux(x_tilde5,gamma1);

plot(x_tilde,y_exp_dim,'o',x_tilde5,y_calc_plot) 

xlabel({'$$\tilde{X}$$'},'Interpreter','latex')

ylabel({'$$\tilde{J}$$'},'Interpreter','latex')

axis([0 1 0 3]);

legend(['beta=' num2str(beta1)])

end

function flux = getflux (x_tilde,gamma)

flux=[1,1,1,1]; %initial guess

    for i=1:length(x_tilde)

    F=@(x) (exp(x)-(1+x+x.*x./2))./(1+gamma.*exp(-x)); %integral of the function


                fun=@(z) x_tilde(i) - (quad(@(x)(exp(x)-x-1),0,z)-quadgk(F,0,z)); %Formula for x_tilde

                z=fzero(fun,1);

                flux(i)=1/(1+gamma*exp(-z)); 
    end
end            

Now, I want to fit another parameter, z0, which should appear in the fifth line from the end instead of the two zeros as follows:

function [K_dep_fit,z0_fit]=twovar
options = optimset('PlotFcns',@optimplotfval); %Plotting the iteration Process\
[K_dep_fit,z0_fit] = fminsearch(@criterion, [0.0964,1], options); %initial value for K_dep=0.0964 and z0=1 
   function crit = criterion (K_dep,z0_fit)
     
       %%%%%constants%%%%%   
y_exp=[7.595131579, 8.238253968,8.236679198,8.582412281,9.562149123]; 
x_mlad1=1.2e-2;x_mlad2=1.7e-2;x_mlad3=2.2e-2;x_mlad4=2.7e-2;x_mlad5=3.2e-2;x_mlad=[x_mlad1, x_mlad2, x_mlad3, x_mlad4, x_mlad5]; 
V_exp=100; 
V=V_exp/60*10e-9; 
area = 289632*1e-12; OD =0.02; C0 = OD * 2.56e14  ;a2 = 0.52e-6; mu=8.9e-4;
g=9.81; deltap=51.793;Vs=2/9*deltap*g*a2*a2/mu; Vs2=6.003e-8; w=0.5e-2;
h=0.03e-2;l=5.3e-2;sigma=6*V/(w*h*h); sigma2=180.84;Kb=1.38064852e-23; T=293.15; 
D=Kb*T/(6*pi*a2*mu); D2=3.897e-13; 

k1=K_dep*D/a2; 
beta1=Vs./k1; 
gamma1=beta1-1; 
x_tilde=x_mlad.*Vs.^3./(sigma.*D.^2);
y_exp_dim=y_exp./area./80./C0./Vs;  

    y_calc=getflux(x_tilde,gamma1); 
    crit=sum(((y_exp_dim-y_calc)./y_exp_dim).^2);    

end

x_tilde5=0.00001:0.01:4;

y_calc_plot=getflux(x_tilde5,gamma1);

plot(x_tilde,y_exp_dim,'o',x_tilde5,y_calc_plot) 

xlabel({'$$\tilde{X}$$'},'Interpreter','latex')

ylabel({'$$\tilde{J}$$'},'Interpreter','latex')

axis([0 1 0 3]);

legend(['beta=' num2str(beta1)])

end

function flux = getflux (x_tilde,gamma)

flux=[1,1,1,1]; %initial guess

    for i=1:length(x_tilde)

    F=@(x) (exp(x)-(1+x+x.*x./2))./(1+gamma.*exp(-x)); %integral of the function


                fun=@(z) x_tilde(i) - (quad(@(x)(exp(x)-x-1),z0,z)-quadgk(F,z0,z)); %Formula for x_tilde

                z=fzero(fun,1);

                flux(i)=1/(1+gamma*exp(-z)); 
    end
end              

Basically, K_dep appears in the function ''crit''. but z0 appears only in ''getflux'' function. I've attached a flow chart as a JPEG for the one parameter fitting and for the two parameters fitting (as far as I undertand it).

The new code obviously doesn't work beacuse there must be some changes with array dimensions due to seeking for 2 parameters and not just one.

I would be very happy to get some help from you guys.

Thank you,

Roi