F = 1./sqrt(f)+4*log((E./3.7)+(1.256./Re*sqrt(f)));
Your initial conditions are scalar, so your f will be scalar, so 1./sqrt(f) will be scalar
Your E is a vector of length 7, so 4*log((E./3.7) will be a vector of length 7.
Your Re is a vector of length 9999000, so (1.256./Re*sqrt(f))) will be a vector of length 9999000 .
You cannot add a vector of length 7 to a vector of length 9999000.
If you are trying to solve over every combination of Re and E then you will need to do something like:
[gE, gRe] = ndgrid(E, Re);
F = arrayfun(@(G, RE) fsolve(@(f) colebrook(f, G, RE), 0.1), gE, gRe);
Unfortunately, the function is nonlinear, and needs to be solved individually.
You should consider, though, that the solution increases with Re and decreases with E, so you might be more efficient looping, moving to increasing value of one parameter and using the previous result as a guess.
