[B,spc] = eprload('4_5mW 9K.DTA');
B = B/10;
plot(B,spc);
% Set up experiment
Exp.mwFreq = 9.64;
Exp.Range = [min(B) max(B)]; % mT
Exp.nPoints = numel(spc);
% Set up spin system with starting parameters
Sys0.g = [1.94 4.65 5.75];
Sys0.gStrain = [0.0351 0.9876 0.6111];
% Set up Vary structure with all parameters that can be varied
SysVary.g = [1 1 1]*1.9;
SysVary.gStrain = [1 1 1]*1.9;
% Function handles: needed for esfit()
f = @sin;
f(0.5)
g = @(x)x^2;
g(5)
% Call esfit
esfit(@pepper,spc,Sys0,SysVary,Exp);
I would refer to https://doi.org/10.1016/S0022-2364(78)80015-8 for the method of extracting E/D from a rhombic g-tensor fit, but a warning, that method is only really valid when |D| is much greater than the Zeeman interaction.
Alternatively you could fit your spectrum using its actual spin multiplicity and zerofield splitting and extract E/D directly that way.