As suggested in a previous answer, numerical integration is a solution in this case.
Below you will find 2 ways of achieving this using R.
The first solution uses an approximation of your kernel density estimates using spline interpolation. A second alternative is to use a trapezoid approximation (one easy to use function can be find in the R pkg sfsmisc).
The term $f(x)$ is computed for the "x"-values in the est object i.e. using dnorm(est$x)
Hope this helps.
set.seed(1)
x1 = rnorm(250)
#Computing of bandwidth
h1 = bw.ucv(x1)
#Estimate
est = density(x1,kernel = "gaussian",bw = h1)
# Create approx func obj and integrate()
splxy = splinefun(est$x, (est$y - dnorm(est$x))^2)
integrate(splxy, lower = min(x1), upper = max(x1)) # res 0.0006681713
# Simple numerical integration in pkg sfsminsc
sfsmisc::integrate.xy(x = est$x, (est$y - dnorm(est$x))^2) # res 0.0006682952