Your ep is 1 x 11 so for each scalar x your f would want to return 1 x 11. integral() does not expect that: integral() expects that for a scalar x, the function would return a scalar. You can use the 'arrayvalued' option to integral to permit a vector to be returned.
However if you do that, then note you overwrite all of P each iteration of the for loop. You could switch to assigning to P(i) but then you would have the problem that that would be a scalar location. So you would need to assign to P(i,:) . But then you have
pr(i) = - (3/2) .* mu .*q .* P
which is outside the loop and would try to assign P to the scalar location pr(i) where i has the last value it was assigned, which would be length(ep) . With P being a 2D array, assigning to a scalar location pr(i) is not going to do well. And then you should have another look at what you are plot()'ing.
It is an odd coincidence that you loop to length(ep) and that you are producing a vector of length(ep) in each iteration.
I would suggest that at the very least you should be changing
f = @(x) 1./(ho*(1 + (ep .* sin(2 * pi * x /L)))).^3;
to
f = @(x) 1./(ho*(1 + (ep(i) .* sin(2 * pi * x /L)))).^3;
If you do that then a bunch of what I said does not apply any more... but the part about overwriting all of P on each iteration does.
