differential equations - Numerically integrating solution obtained from NDSolve method

In the following example, $u(x)$ is found numerically using NDSolve method.

 F = 1/1000

    h = 12000/1000
    d = 10/10
    L = 1000
    W = 3
    phi[x_] := 
     Piecewise[{{(1/2)*(1 - Tanh[((L*x)/(d))]), 
        x <= 1/2}, {(1/2)*(1 + Tanh[((L*(x - L/L))/(d))]), x > 1/2}}]
    vE[x_] := x*(1 - x)*4
    s = NDSolve[{u''[x] == (h*L*L/(d*d))*phi[x]*phi[x]*u[x] - 
         F*L*L*(1 - phi[x]), u[-W*d/L] == 0, u[1 + W*d/L] == 0}, 

      u, {x, -W*d/L, 1 + W*d/L}, Method -> "StiffnessSwitching", 
      WorkingPrecision -> 40, InterpolationOrder -> All]
    diff[x_] := (u[x] - vE[x])*(u[x] - vE[x])
    Plot[Evaluate[{diff[x]} /. s], {x, W*d/L, 1 - W*d/L}, 
     PlotRange -> All]

Which works perfectly. I need to see what is mean square error between obtained solution and another function $vE(x)$.

sum = 0;
Do[
 first = W*d/L;

 second = 1 - W*d/L;
 {sum = sum + diff[first + (i/100)*(second - first)]},
 {i, 0, 100, 1}]
Evaluate[sum]

but this gives only expression but not value. I think this is because $u(x)$ is obtained at discrete points only and is not defined on the points on which I have calculated error. I also tried using integration,

intVal = NIntegrate[({u[x]} /. s - vE[x])*({u[x]} /. s - vE[x]), {x, 
   W*d/L, 1 - W*d/L}]

but this gives long error message ending with,

 "...is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing"

How can I evaluate this integral?

Answer

You can integrate the mean square error mse at the same time as computing u[x].

s = NDSolve[{
    u''[x] == (h*L*L/(d*d))*phi[x]*phi[x]*u[x] - F*L*L*(1 - phi[x]),
    u[-W*d/L] == 0, u[1 + W*d/L] == 0,
    mse'[x] == (u[x] - vE[x])^2, mse[-W*d/L] == 0},
   {u, mse}, {x, -W*d/L, 1 + W*d/L}, Method -> "StiffnessSwitching", 

   WorkingPrecision -> 40, InterpolationOrder -> All];

 Plot[Evaluate[mse[x] /. s], {x, W*d/L, 1 - W*d/L}, PlotRange -> All]

You seem to be interested in this change:

mse[1 - W*d/L] - mse[W*d/L] /. First@s

(* 8198.070964448656291179158359833465311096 *)

---

The problem with the code

 intVal = NIntegrate[({u[x]} /. s - vE[x])*({u[x]} /. s - vE[x]), {x, W*d/L, 1 - W*d/L}]

is a syntax issue. The expression that /. tries to apply is s - vE[x]. This can be seen from the TreeForm of the expression:

Hold[({u[x]} /. s - vE[x])] // TreeForm

In other words, ReplaceAll has lower precedence than Plus (represented by the minus sign). The proper code is

intVal = NIntegrate[((u[x] /. First@s) - vE[x])^2, {x, W*d/L, 1 - W*d/L},

  WorkingPrecision -> 40]

(* 8198.070964448656291179224357022321689725 *)