A couple of possibilities for reducing the execution time:
1) Use the 'roots' function instead of 'fzero' and select the single root of the five whose imaginary part is zero for this polynomial.
2) Do a sort on u and use the y solution for one sorted u value as an initial estimate of y for the next higher u. Because successive u values will necessarily be close together, these estimates should already be fairly accurate. Also it might help to use the Newton-Raphson method instead of 'fzero' in combination with this - you already have an expression for the derivative of F(x), namely f(x). Since you are only doing a histogram of y it shouldn't matter that these are also now in sorted order, but if necessary their original order could easily be restored.
Note: The 'find' operation is unnecessary here since all y values will lie between 0 and 1, which is true because F(0)=0 and F(1)=1 and F is monotone increasing in between.
