Well, I'm not saying that this is the best way, but you could try this:
n = length(fInf);
m = length(fLin);
fLin_repeated = repmat(fLin, [1, n]);
S2_repeated = repmat(S2, [1, n]);
fInf_repeated = repmat(fInf', [m, 1]);
fSup_repeated = repmat(fSup', [m, 1]);
tmp1 = (fLin_repeated > fInf_repeated) & (fLin_repeated < fSup_repeated);
P = sum( S2_repeated .* tmp1 )';
I particularly don't like the last line. I feel like there is a function to do the masking part instead of doing the multiplications, but for the life of me I can't remember what it is. The downside is that you will have 4 different 50000 by 1000 matrices, which may or may not be problematic.
Hope this helps.
P.S. Testing it on my machine, I run out of memory on making the 4th matrix. To get around this you could either cut up S2 into smaller chunks (maybe to 10000 at a time instead of all 50000) or you could calculate the tmp1 variable, clear fLin_repeated, fSup_repeated, and fLin_repeated, then generate S2_repeated and P.