Hello guys!
I need your help to find the parameters by fitting the data using a differential evolution algorithm on a non-linear mathematical model containing double integral which can not be determined analytically but numerically. I used multiple functions to reduce the expression of integration,

F = {{\[Lambda]1, 0, 0}, {0, \[Lambda]1, 0}, {0, 0, 1/\[Lambda]1^2}};
c = Transpose[F] . F;
EGreen = 1/2 (c - IdentityMatrix[3]);
bleft = F . Transpose[F];
Nref = {Sin[\[CapitalTheta]] Cos[\[CapitalPhi]], 
   Sin[\[CapitalTheta]] Sin[\[CapitalPhi]], Cos[\[CapitalTheta]]};
ncur = F . Nref;
w1 = D[c1/(2 c2) (Exp[c2 (Imat)^3] - 1), Imat]; 
Imat = Simplify[
   Tr[EGreen. (TensorProduct[Nref, Nref])] + Sqrt[
    Tr[EGreen^2 . (TensorProduct[Nref, Nref])]]];
Hmod = 1/
    2 (TensorProduct[Nref, Nref] + 
      1/(2 Sqrt[
        trEEH]) (EGreen . (TensorProduct[Nref, 
            Nref]) + (TensorProduct[Nref, Nref]) . EGreen)) /. 
   trEEH -> Tr[EGreen . EGreen . (TensorProduct[Nref, Nref])];
hmod = Simplify[F . Hmod . Transpose[F]];

Sigma11BiaxialMatchedInvariant[n_, b_, c1_, c2_] := 
  Table[0.3 (\[Lambda]1^2 - 1/\[Lambda]1^4) + 
    NIntegrate[
     4 n Sqrt[b/(2 Pi)]
        Exp[2 b (Cos[\[Alpha] - \[CapitalPhi]] \
Sin[\[CapitalTheta]])^2]/
       Erfi[Sqrt[2 b]]*2 w1*(-hmod[[3, 3]] + hmod[[1, 1]])*
       Sin[\[CapitalTheta]] /. {\[Alpha] -> 0}, {\[CapitalPhi], 0, 
      2 \[Pi]}, {\[CapitalTheta], 0, \[Pi]}], {\[Lambda]1, 0.5, 1.5, 
    0.04}];

fit = NonlinearModelFit[data, 
  Sigma11BiaxialMatchedInvariantSigma11[n, b, c1, c2], {n, b, c1, c2}, \[Lambda]1, 
  Method -> {NMinimize, 
    Method -> {"DifferentialEvolution", "ScalingFactor" -> 0.9, 
      "CrossProbability" -> 0.1, 
      "PostProcess" -> {FindMinimum, Method -> "QuasiNewton"}}}]

I have attached the data. Please help me!

Attachments:

GST33.mat