Best Curvefit method for Viscoelastic Material

Hello, I currently have a curvefit code in place to curvefit and find the required parameters for experimental data.

Please find attached the data file if needed.

The Code I am using is as following. While I am satisfied with the results, Can anyone help me in knowing if this is the best possible way to do it.

See attached code.

if true
  if ~exist('Processed','dir')
          mkdir('Processed')
end
Files = dir('*.xls') 
numfiles = length(Files);
for num = 1:numfiles 
Data = xlsread(Files(num).name);
time=Data(:,1);
load_data=Data(:,5);
depth=Data(:,2);
[M,I] = max(load_data(:));
tR=time(I)+0.01;
time=time(I:length(time));
contact_depth=depth(I);
load_data=load_data(I:length(load_data));
lowest_load = load_data(end);
k=0.0025;
R=2.5e-06;
x0=[1e6 4e6 1e8];
%  x0=[2e4 4e4 1e6]; % Initial design point
xl=[-inf -inf -inf];
xu=[inf inf inf];
sample=strtok(Files(num).name,'.');
[fitresult,gof,fig]=createFit(time,load_data,tR,k,R,x0,sample,xl,xu);

E1out=fitresult.E1; E2out=fitresult.E2; etaout=fitresult.eta; R2=gof.rsquare; xls_data{1,1}='E1= '; xls_data{2,1}='E2= '; xls_data{3,1}='eta= '; xls_data{4,1}='R2= '; xls_data{5,1}='Peak Load= '; xls_data{6,1}='Einf= '; xls_data{7,1}='Z= '; xls_data{8,1}='Lowest load='; xls_data{1,2}=E1out; xls_data{2,2}=E2out; xls_data{3,2}=etaout; xls_data{4,2}=R2; xls_data{5,2}=M; Einf=((E1out)*(E2out))/((E1out)+(E2out)); xls_data{6,2}=Einf; xls_data{7,2}=contact_depth; xls_data{8,2}=lowest_load;

xlswrite([strcat('\Processed\',Files(num).name)], xls_data, 1, 'A1'); % xlswrite([strcat('\\Output\Optimized_Values_',sample ,'.xls')], xls_data, 1, 'a1') end end