I want to fit the data with the mathematical model containing numerical integration, please help me!

I have the data as GST33.mat file. I want to fit this data with the mathematical model where ; and to get parameters . I tried very hard to get parameters but I'm getting errors. Please help me.

This is the code which I tried.

%AIfit.m
function par = AIFit
    load('GST33.mat','Expression1');
    xdata = Expression1(:,1);
    ydata = Expression1(:,2);
     % Inital guess for parameters:
    n0= 1;
    c10=1;
    c20= 0;
    b0=1;
    par0 = [n0;c10;c20;b0];
     % lsqcurvefit is in the optimization toolbox. 
     % fit, from the curve fitting toolbox may be an alternative
     par = lsqcurvefit(@Kdvec,par0,xdata,ydata);
     % Check fit:
     figure;
     plot(xdata,ydata,'or')
     hold on
     xplot = linspace(min(xdata),max(xdata));
     plot(xplot,AIdvec(par,xplot))
     hold off
end
%Kdvec.m    
function AIvec = Kdvec(par,xdata)
  AIvec = zeros(size(xdata));
for i = 1:length(xdata)
    AIvec(i) = Kd(par,xdata(i));
end
%Kd.m
end
 function K = Kd(par,lambda3)
    lambda = @(theta,lambda3) sqrt(((lambda3)^2)*((cos(theta)).^2)+((lambda3)^(-1)) ...
                     *(sin(theta)).^2);
    rho    = @(theta,lambda3,par) 4*sqrt(par(4)/(2*pi))*(exp(2*par(4)*cos(theta).^2)./erfi(sqrt(2*par(4))));
    Dw     = @(theta,lambda3,par) 2*par(2)*lambda(theta).*(-1+lambda(theta).^2).*(exp(par(3)* ...
              (lambda(theta).^2-1)).^2);
    integrand = @(par,lambda3, theta) pi*par(1)*rho(theta,lambda3,par).*Dw(theta,lambda3,par).*(2*(lambda3^2)*(cos(theta).^2)-(lambda3^(-1)) ...
                .*(sin(theta).^2))./lambda(theta,lambda3)+0.3*(lambda3^2-lambda3^(-1));
    K =  integral(@(theta) integrand(par,lambda3,theta),0,pi);
end   

Answer 1

[Torsten]

%AIfit.m
%function par = AIFit
    load('GST33.mat','Expression1');
    xdata = Expression1(:,1);
    ydata = Expression1(:,2);
     % Inital guess for parameters:
    n0= 1;
    c10=1;
    c20= 0;
    b0=1;
    par0 = [n0;c10;c20;b0];
     % lsqcurvefit is in the optimization toolbox. 
     % fit, from the curve fitting toolbox may be an alternative
     par = lsqcurvefit(@Kdvec,par0,xdata,ydata);
     % Check fit:
     figure;
     plot(xdata,ydata,'or')
     hold on
     xplot = linspace(min(xdata),max(xdata));
     plot(xplot,AIdvec(par,xplot))
     hold off
%end
%Kdvec.m    
function AIvec = Kdvec(par,xdata)
  AIvec = zeros(size(xdata));
for i = 1:length(xdata)
    AIvec(i) = Kd(par,xdata(i));
end
%Kd.m
end
 function K = Kd(par,lambda3)
    lambda = @(theta,lambda3) sqrt(((lambda3).^2)*((cos(theta)).^2)+((lambda3).^(-1)) ...
                     .*(sin(theta)).^2)
    rho    = @(theta,lambda3,par) 4*sqrt(par(4)/(2*pi)).*(exp(2*par(4).*cos(theta).^2)./erfi(sqrt(2*par(4))))
    Dw     = @(theta,lambda3,par) 2*par(2).*lambda(theta,lambda3).*(-1+lambda(theta,lambda3).^2).*(exp(par(3)* ...
              (lambda(theta,lambda3).^2-1)).^2)
    integrand = @(par,lambda3,theta) pi*par(1).*rho(theta,lambda3,par).*Dw(theta,lambda3,par).*(2*(lambda3^2).*(cos(theta).^2)-(lambda3^(-1)) ...
                .*(sin(theta).^2))./lambda(theta,lambda3)+0.3*(lambda3^2-lambda3^(-1))
    K =  integral(@(theta) integrand(par,lambda3,theta),0,pi);
end