Fitting Inverse Power Polynomial to dataset with very small values

I'm trying to fit an inverse power polynomial of 4th power, but the fit function is not fiiting the data correctly.

Even when I manually give the correct parameter estimates, it fails to reproduce the fit. (Just finding the parameters manually is not enough, as I need a generic code which will work with different datasets).

I've been trying to play around with the fit options to fix the fit, but I didn't succeed.

It seems like the issue is that I am dealing with very small numbers (1e-18) and coefficients which are orders of magnitude different, and Matlab is having trouble to fit it.

I would appreciate any kind of advice, as I've been working on this for quite long and didn't get any meaningful fittings.

I attached a dataset, and copied a functioning example below.

fid = fopen('phase05.txt', 'r'); % pointer to file
data = []; % initialize; don't know size of matrix;
row = 1; % initialize;
while ~feof(fid) % loop until file ends;
    aline = fgetl(fid); % read in a new line;
    data(row,:) = sscanf(aline, '%g');  %convert to number;
    row = row + 1;
end
fclose(fid); % close file;

f = data(:,1); % extracting the values of frequency
noise = data(:,2); % extracting the values of logarithmic power spectral density of frequency L(f)

Sphi = 2 * (10.^(noise./10)); % converting from L(f) to Sphi(f) 
fitPLaw = @(f) 2*p0 + (2*p1)./f + (2*p2)./(f.^2) + (2*p3)./(f.^3) + (2*p4)./(f.^4);

myFitFunc = fittype('2*p0 + (2*p1)./f + (2*p2)./(f.^2) + (2*p3)./(f.^3) + (2*p4)./(f.^4)', ...
    'dependent','Sphi','independent','f', ...
    'coefficients',{'p4', 'p3', 'p2', 'p1', 'p0'});

[myFit, goffit, outputF] = fit(f, Sphi, myFitFunc, 'StartPoint', [10^(-18),10^(-18),10^(-18),10^(-18),10^(-18)],...
    'DiffMinChange',1e-23,'DiffMaxChange',1e-4,'TolFun',1e-50,'MaxIter',10000000,'Robust','Bisquare',...
    'TolX',1e-20,'MaxFunEvals',10000);

pCoeff = fit(f, Sphi, myFit,'MaxIter',1000); %, 'StartPoint', [1e-8, 1e-8, 1e-8, 1e-8, 1e-8], 'Method', 'LinearLeastSquares');
pArray=coeffvalues(pCoeff); 

noiseSim=2*pArray(5) + (2*pArray(4))./f + (2*pArray(3))./(f.^2) + (2*pArray(2))./(f.^3) + (2*pArray(1))./(f.^4);
myNoise=2*1e-18 + (6*1e-12)./f + (2*1e-11)./(f.^2) + (2*2e-5)./(f.^3) + (2*1e-2)./(f.^4); % phase05
logF=log10(f);          % log of frequency
logPhi=log10(Sphi);     % log of the data
logFit=log10(noiseSim); % log of the extracted fit on Sphi

figure;
plot(logF,logFit,logF,logPhi,logF,log10(myNoise));
legend('L-fitted,log(f)','L-data,log(f)','using values manually ,log(f)')