You can do something like this:
y_fcn = @(CT,x) CT(1) * 74.826 * (x./CT(4)).^3 * integral(t.^4.*exp(t)./(exp(t)-1).^2), 0, CT(4)./x) + ...
CT(2) * 24.942 * (CT(5)./x).^2 .* exp(CT(5)./x)./(exp(CT(5)./x)-1).^2 + CT(3)*24.942*(CT(6)./x).^2.*exp(CT(6)./x)./(exp(CT(6)./x)-1).^2;
err_fcn = @(CT,x,y) sum((y-y_fcn(CT,x)).^2);
CT0 = [1 2 3 4 5 7]; % Some sensible initial guess for C and T, that we concatenate together as CT = [C,T]
CT_best = fminsearch(@(CT) err_fcn(CT,x_obs,y_obs),CT0)
You most likely have to use element-wise operations in the equation for y (I might have missed one or two multiplictions or power-operations). Then it should be fine to use ordinary least-square fitting of the C and T parameters (provided you have enough data-points, and they have the same standard deviation - but I expect you to know this already).
HTH
