The goal is to minimise a functional in this case I want to find the input function, defined by its (many) values on a grid for which the functional attains its minimum. Below is a piece of code where I i) make a grid, ii) evaluate some function on the grid (that will serve as the initial guess) iii) define some functions to interpolate and integrate iv) try to minimise the objective function (called total energy). All user defined functions seem to work fine, but when I evaluate FindMinimum NINtegrate throws error:

Integrate::inumr: ... has been evaluated to non-real values.

I read somewhere (on stack) that I should define an array to let Mathematica know I only want it to look for Real values on the grid points. I looked at this related question: http://mathematica.stackexchange.com/questions/4002/minimizing-a-function-of-many-coordinates. (Note: I have limited the grid to 5 points, but I will in fact need 40000 to get accurate integrated values.)

   genrhoi = Function[{rho1, h}, N[Range[rho1, 0.3, h]]]; 
    grid = genrhoi[0., 1./16.] + 10^(-6)
    densityongrid = (Exp[-3.*Abs[#1]] & ) /@ grid
    setupintegrand = Function[{functionongrid}, functionongrid*grid^2.]; 
    Integrator = Function[interpolatedintegrand, 4.*Pi*NIntegrate[interpolatedintegrand[x], 
          {x, grid[[1]], grid[[-1]]}, WorkingPrecision -> 10]]; 
    Interpolator = Function[interpolatethis, Interpolation[
         Transpose[{grid, interpolatethis}], InterpolationOrder -> 4]]; 
    nuclearpotential = Function[{position, charge, grid}, 
        (-1.*(charge/EuclideanDistance[position, #1]) & ) /@ grid]; 
    externalpotential = nuclearpotential[0., 7., grid]; 
    VRHO = Function[density, Integrator[Interpolator[setupintegrand[
           Evaluate[externalpotential*density]]]]]; 
    TF = Function[density, Integrator[Interpolator[setupintegrand[
           Evaluate[density^(5./3.)]]]]]; 
    totalenergy := TF[#1] + VRHO[#1] & ; 
    vars = Array[a, Length[grid], 1]
    FindMinimum[totalenergy[vars], Transpose[{vars, densityongrid}], MaxIterations -> 5, 
      PrecisionGoal -> 4, Gradient -> Evaluate[externalpotential + vars^(2./3.)], 
      Method -> "QuasiNewton", EvaluationMonitor -> Print[totalenergy[vars], vars, 
        ListLinePlot[Transpose[{grid, vars}]]]]