By the way, when I copy the code from the browser and paste it into Mathematica (I am nostalgic of the old Greek name), the integral from 1 to Infinity becomes from Infinity to 1. It looks like a bug.

Define epsilon to be a function of r (ε[r_] :=...) and change the formula for f[r] to use ε[r].

When the function f[r1, P1] is evaluated, first the literal instances of the symbols r and P in the (unevaluated) formula for f are replaced by r1 and P1. The symbol epsilon remains the symbol epsilon. After this replacement, the application of f is complete, and the resulting expression is evaluated. During this evaluation, epsilon is evaluated and the symbol r in it appears. But the application of f is done; so the r remains, the value of the function not a number.

A standard check one should try when a plot is blank is to evaluate the function at some numbers. For instance fPlot[0.1]. If you don't get a number, then that's why the plot is blank. The next step is to track down why the value was not a number.