Hi Zeyad,
Here is an alternative function quad() with its set tolerance to compute integration that changes the solutution a bit.
clearvars
length_shaft = 0.4;
% Determining 4 values of beta
f = 'cos(x)*cosh(x)-1';
OPTs = optimset('tolX',1e-17);
Roots = [fzero(f,4.0, OPTs) fzero(f,7.0, OPTs) fzero(f,10.0, OPTs) fzero(f,14.0, OPTs)];
beta = Roots./length_shaft;
% Determining Constant (T) (a value for each beta)
for ii=1:numel(beta)
B(ii) = quad(@(y) ((cos(beta(ii)*y) + cosh(beta(ii)*y)) - (((cos(beta(ii)*length_shaft) - ...
cosh(beta(ii)*length_shaft))/(sin(beta(ii)*length_shaft) - ...
sinh(beta(ii)*length_shaft)))*(sin(beta(ii)*y) + sinh(beta(ii)*y)))).^2, 0, length_shaft, 1e-16);
T(ii) = (B(ii)*2.5)^(-0.5);
end
format long
B = double(B)
T = double(T)
An alternative way would be for the integration part to write your own code based on Simpson's rule that gives much higher accuracy.
Good luck.
