Every multiplication in ‘Fermi’ needs to be vectorised (with element-wise operations), the Symbolic Math Toolbox is not necessary, (so replace solve with fsolve), and use the ArrayValued name-value-pair in the integral call.
With those changes, your code:
dpsi = -5;
etaK = 2;
etaL = 4;
j0 = 1;
UT = 1;
Fermi = @(eta) sqrt(pi) ./ (2 .* (3/4 .* sqrt(pi) .*( eta.^4 + 33.6.*eta .* (1.0-0.68.*exp(-0.17*(eta+1).^2) ) + 50).^(-3/8) + exp(-eta) ));
InnerInt = @(eta, flux) 1/ ( flux/Fermi(eta) + dpsi/UT);
integralEq = @(flux) integral( @(eta) InnerInt(eta, flux), etaK, etaL, 'ArrayValued',1);
jKG = j0 * fsolve(@(x) integralEq(x) - 1, 10)
now produces this result:
jKG =
-3033.2578125
