@Saeed I don't see any inaccuracy in getting the points from the graph. After all, how is the curve represented in the computer? It is an ordered collection of (x,y) pairs.

@Saeed Get points from the graph. If the graph exists in any computerized form, that should be fairly easy.

@Carl Love My Answer above has been significantly updated.

Nobody wants to retype your integrals. So please post your work as a worksheet. You can use the green uparrow, which is the last item on the second row of the toolbar in the MaplePrimes editor.

@Chris_45 

I have tried Statistics:-NonlinearFit on many different interpolating models, including sines, exponentials, rational functions, and products and compositions of these. I can't get any reasonable fit. The "associated function" for my ArrayInterpolation curve is a piecewise cubic spline with very many pieces. I doubt that that would be of much use to you, but if you want it, just use Preben's Spline command.

Part of the problem is that your independent variable, y_data, is very unevenly spaced: It has huge gaps in the middle and is very tightly spaced at the high end.

This is just an aside, not an answer to your Question. If you know that your sums can't be evaluated, you should use Sum instead of sum. If you use the latter, the system wastes time trying to evaluate and re-evaluate the sums everytime the expression containing them is used. Also, one of your transformations may lead to an unintended simplification.

@brian bovril 

Could you provide a more-standard reference to the paper so that I can find it?

@Carl Love 

I don't know much about Maplels, but I found the following in the help page ?Maplets,Elements,TextField: The command Maplets:-Elements:-TextField has an option tooltip= "...".

@Alejandro Jakubi Thanks for spotting. Code corrected above.

@JohnS Pattern matching "any name except x" is possible like this:

define(
     L, orderless, multilinear,
     L(a::And(name,Not(identical(x)))*b::function, c::function) = a*L(b,c),
     L(
          a::And(name,Not(identical(x)))*b::function,
          c::And(name,Not(identical(x)))*d::function
     ) = a*c*L(b,d)
);

@JohnS Under your definitions for L, we have L(x*f(x)) = x*L(f(x)). If that makes sense for your application, fine; but it doesn't fit my usual conception of linearity.

@JohnS 

It will work if you first put them on the sequence of constants:

constants:= constants, c, d, e, f;

@acer 

It wasn't the extra function calls that worried me; indeed, my procedure above has only one call per dimension, which is trivial. But I was worried about the one-element-at-a-time list building; each tuple is created that way. There might be a recursive way using Vectors. For a Cartesian product we know the number of tuples and their size from the start.

@Chia Go to help page ?parameter_classes and read the ninth section "End-of-Parameters marker, $". This usage of $ has no relation whatsoever to the more-familiar usage of $ as a sequence-building operator.

What do you mean by "high powers, like 10^17 by 10^12"?

Your story seems strange: Raising a 250x250 Matrix of floats to the 10^17 power is nearly instantaneous. What exactly are you doing? Why is a loop needed to raise a Matrix to a power?