I think what you might be missing is that FIND is returning an index, so the eta value you seek is eta(F12_p). But I didn't really try your approach. Here's how I might have done it.
etaspan = -500:0.001:500;
f = @(x) (x.^(1/2))./(1+exp(x-etaspan)); % A vector-valued function.
y = integral(f,0,500,'ArrayValued',true); % Calculate the integrals in one call.
% If y is already sorted, just use ys = y and idx = 1:length(etaspan).
% Note that we do assume that the y values are unique.
[ys,idx] = sort(y);
% Use the table to look up the nearest eta value given a y value.
findeta = @(y)etaspan(interp1(ys,idx,y,'nearest','extrap'));
FIND uses a linear search. INTERP1 uses a binary search, which should be faster. Now if you didn't want to use the table, fzero might work.
f = @(x,eta)(x.^(1/2))./(1+exp(x - eta));
findeta = @(y)fzero(@(eta)integral(@(x)f(x,eta),0,500) - y,0);
