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?