How to make a function depend on the solution of an integration?

I have two functions which are the linear and non-linear solution of an equation. I'm trying to calculate where does the linear solution is a good solution of my equation, when compared with the non-linear. To do so I am using

$$\int_{\tfrac{a}{2}}^\infty\left(\sigma^{NL}(x)-\sigma^{L}(x)\right)^2\mathrm dx$$

The variables I'm using are

T = 298.15(*K*);
k1 = 1.38064852*10^-23 (*J/K*);
ϵ0 = 8.85418781761*10^-12 (*C^2/N m^2||F/m||C/V m*);
ϵ = 78.5;
ϵR = ϵ ϵ0 (*C^2/N m^2||F/m||C/Vm*);
c = 2.99792458*10^8(*m/s*);
z = 1;
e1 = 1.60217733*10^-19 (*C*);
NA = 6.022140857*10^23 (*mol^-1*);
a = 4.25(*Å*);
a1 = a*10^-10 (*m*);
x1 = x*10^-10 (*m*);

And the main code is

nMKS[ρ0_] = ρ0 *10^3 NA (*m^-3*);
κMKS[ρ0_] = Sqrt[(2 nMKS[ρ0] (e1^2) (z^2) )/(ϵR k1 T)](*m^-1*);
LσMMKS[ρ0_, ψH_] = ϵR κMKS[ρ0] ψH/1000  (*C/m^2*);
LψPBMKS[ρ0_, ψH_, x_] = ψH Exp[κMKS[ρ0] (a1/2 - x1)] (*mV*);
LσMKS[ρ0_, ψH_, x_] =
Piecewise[{
{0, 0 <= x1 < a1/2},
{ϵR κMKS[ρ0] LψPBMKS[ρ0, ψH, x]/1000, x1 >= a1/2}
}](*C/m^2*);
NLσMMKS[ρ0_, ψH_] = (2 ϵR κMKS[ρ0] k1 T)/( z e1) Sinh[(z e1 )/(2 k1 T) ψH/1000] (*C/m^2*);
ZMKS[ψH_] = (z e1 ψH/1000)/(2 k1 T);
NLσMKS[ρ0_, ψH_, x_] =
Piecewise[{
{0, 0 <= x1 < a1/2},
{(4 ϵR κMKS[ρ0] k1 T)/(z e1) (((Exp[2 ZMKS[ψH]] - 1) Exp[κMKS[ρ0] (a1/2 - x1)])/(
(Exp[ZMKS[ψH]] + 1)^2 - ((Exp[ZMKS[ψH]] - 1) Exp[κMKS[ρ0] (a1/2 - x1)])^2)), x1 >= a1/2}
}](*C/m^2*);

Hence, when I try to integrate my equation with

Dif2σ[ρ0_, ψH_, x_] = FullSimplify[(NLσMKS[ρ0, ψH, x] - LσMKS[ρ0, ψH, x])^2]
Vσ[ρ0_, ψH_] = FullSimplify[Sqrt[Integrate[Dif2σ[ρ0, ψH, x], {x, a/2, Infinity}]]]

Mathematica is not able to solve it. How can I force Mathematica to actually solve this and thereafter define a function from its solution?