The result of interpolating a table is not numeric

Question

I'm having trouble in using interpolating functions with Mathematica. Some interpolating functions I defined, say $f[x]$, if evaluated in a specific point $x_0$ correctly return the numeric value $N[f[x_0]]$, while some others return a non-numeric value $\{f[x_0]\}$, as if the interpolating function were a list of array-like values. This leads me to problems in the subsequent manipulation of these interpolating functions since, for example, when integrating such functions one can get error messages like "Integrand is not numerical at {x} = {...}". I'm now going to show the problem in detail.

I'm using Mathematica to numerically solve a differential equation $d\chi/d\phi=F(\phi)$, with a given function $F(\phi)$, to find $\chi(\phi)$. Then I want to reverse the newly found relationship to obtain $\phi(\chi)$. The problem I'm facing is independent of the complexity of the function $F(\phi)$, so from now on I will just assume $F(\phi)\equiv1$.
Then we are just numerically solving $\,\chi(\phi)=\int d\phi\,$ to find that $\chi(\phi)=\phi$ and reverse this to obtain $\phi(\chi)=\chi\,$.
The code I wrote is the following:

xMax = 40 ;
STEPS = 10000 ;
dx = xMax/STEPS;
CHI[x_] = y[x] /. First[NDSolve[{D[y[x], x] == 1, y[0] == 0}, y, {x, 0, xMax}] ] ;
chi = Interpolation[Table[{{x}, CHI[x]}, {x, 0, xMax, dx}]] ;
phi = Interpolation[Table[{CHI[x], {x}}, {x, 0, xMax, dx}]] ;

Here "CHI[x_]" should be the solution of the differential equation $\,d\chi/d\phi=1\,$. I then defined "chi" to be an interpolating function of "CHI[x_]", namely $\chi(\phi)$. Then I defined the interpolating function "phi" by simply inverting the order of the values I give in the Table command, namely to find $\phi(\chi)$.

Now I would expect "chi" and "phi" to be two interpolating functions, so that by asking for chi[15] and phi[15] in the code I get the numeric values of the two functions in the point "15" (which obviously must both be equal to 15).
By doing this, I get

CHI[15]
chi[15]
phi[15]
Out[7] = 15.
Out[8] = 15.
Out[9] = {15.}

Here is my problem: why both "CHI" and "chi" returns a numeric value 15, while "phi" returns a non-numeric array-like value {15}?

This is critical for my code, because I then have to work with $\phi(\chi)$ to figure out a potential and numerically integrate that potential... and it gives me the error output "Integrand is not numerical at {x} = {...}" for some "{...}" value of x.
In fact, if in this trivial case just shown I ask for NIntegrate the interpolating function "phi", i get the error

NIntegrate[phi[x], {x, 0, chi[xMax]}]
NIntegrate::inum: Integrand InterpolatingFunction[{{0.,40.}},{5,3,0,{10001},{4},0,0,0,0,Automatic,<<3>>},{{0.,0.004,0.008,0.012,0.016,0.02,0.024,0.028,0.032,0.036,<<9991>>}},{{{0.}},{{0.004}},{{0.008}},{{0.012}},{{0.016}},{{0.02}},{{0.024}},{{0.028}},{{0.032}},{{0.036}},<<9991>>},{Automatic}][x] is not numerical at {x} = {0.318293}.
NIntegrate::inum: Integrand InterpolatingFunction[{{0.,40.}},{5,3,0,{10001},{4},0,0,0,0,Automatic,<<3>>},{{0.,0.004,0.008,0.012,0.016,0.02,0.024,0.028,0.032,0.036,<<9991>>}},{{{0.}},{{0.004}},{{0.008}},{{0.012}},{{0.016}},{{0.02}},{{0.024}},{{0.028}},{{0.032}},{{0.036}},<<9991>>},{Automatic}][x] is not numerical at {x} = {0.318293}.
Out[10] = NIntegrate[phi[x], {x, 0, chi[xMax]}]

While the N-integration of "chi" (which has numeric values) correctly returns

NIntegrate[chi[x], {x, 0, xMax}]
Out[11] = 800.

being $\int_0^{40}x\,dx=800\,$.

What is the problem? Why am I not having numeric values in $\phi(\chi)$?