You are asking for normalization. When you do that, coeffvalues are the coefficients for the normalized independent variable, not for the original variable. However, when you build fsvar you are not taking that into account.
It is not exactly that fsvar itself is wrong in itself, but you need extra steps:
meanx = mean(Temp);
stdx = std(Temp);
fsvar = subs(fsvar, x, (x-meanx)/stdx);
Or you could take a different approach that does not use the curve fitting toolbox:
syms x
% reading the data
filename='Results1.xlsx';
sheet='Sheet1';
fs = xlsread(filename,sheet,'A2:A57');
Temp = xlsread(filename,sheet,'B2:B57');
[p, ~, mu] = polyfit(Temp, fs, 4);
fsvar = subs(poly2sym(p, x), x, (x-mu(1))/mu(2));
fss = fsvar^2/(1-fsvar)^3;
ht = matlabFunction(fss);
s1=integral(ht,831,922);
