Fitting custom function involving integrals

Alex Sha anwered a similar question from 'mz86' on May 24th, 2019. I tried to enter my question there but it was taken as my answer to the original question. It was not my answer. My apology.

Here is my concern for which I seek your help.

Alex Shaw's parameter determination for this problem is very different from what 'mz86' obtained from the Matlab curve fitting tool. Why is it so? I too find problems in fitting my custom equations that involve integrals. Also, Alex Shaw only provides the parameters, not the 95% confidence bounds for each of the parameters. What is the difference in the two approaches that yield drastically different parameters for fitting the same data using the same equation? How can I replicate the fits and parameter results obtained by Alex so that I could try to implement the same approach to my curve fitting problem?

Your early response will be gratefully appreciated.

Thank you.

Here is one of my custom equation to fit 1048 data points with errors. Excel data file is attached.

y=VPCexp(x,Om0,299840/H0,a)+5*log10(1+x)+25

function [y] = VPCexp(z,Om0,R0,a)

n = length(z);

y=zeros(size(z));

for i=1:n

zl=z(i);

fun = @(x) (Om0.*((1+x).^3.*exp(a./3.*(1-(1+x).^(-3)))-exp(2.*a./3.*(1-(1+x).^(-3))))+(1+x).^2).^-0.5;

yl = integral(fun,0,zl);

dP=R0*sinh(yl/exp(-a/3*(1-(1+zl)^(-3))));

y(i)= 5*log10(dP);

end

The solution I get is:

General model:

f(x) = VPCexp(x,Om0,299840/H0,a)+5*log10(1+x)+25

Coefficients (with 95% confidence bounds):

H0 = 71.24 (70.47, 72.01)

Om0 = 0.5749 (0.3704, 0.7793)

a = 0.8509 (0.6572, 1.045)

Goodness of fit:

SSE: 1034

R-square: 0.997

Adjusted R-square: 0.997

RMSE: 0.9948