Steven: beautifully simple code on your previous post.  I just fear that MPL is beyond the reach of mortals such as myself because I could never come up with that solution.

If people are in the mood for this kind of stuff, let me post on a set of questions coming up in my libraries:
1) QUESTION 1: MANAGING DERIVATIVES FOR FUNCTIONS
* I am doing a lot of work with functions that may have an associated derivative analytically or numerically determined.
* This is the tip of the iceberg in associating a bunch of information with a mathematical function that all needs to be bundled.
* I have been creating functors like the following:
template <class ParamType>
struct square{
double operator()(double x, const ParamType& params ) { return x*x;}
double gradient(double x, const ParamType& params ){return 2 * x;}
boost::function<double (double,const ParamType& params)>  gradient(){return boost::bind(&square::gradient, *this, _1, _2)}; //I forget the exact notation for bind with member functions, but the idea is that it would return a boost::function.
}; //Analytical derivative

But I would also love the ability to use finitedifference or autodifferentiation algorithms for more complicated examples.  Perhaps something like:

struct square : finite_differences<int Order>{
double operator()(double x, const ParamType& params ) { return x*x;}
}; //Analytical derivative

* Ideally here, finite_differences<int Order> would be able to generate the derivative using order taylor series approx and the operator() to evaluate the function.  Even here I am not sure the best way to organize it.

2) QUESTION 2: ADAPTATION AND USE LIKE BOOST::FUNCTION
Now, lets say that there was something similar to boost::function but which had the "derivative" concept to call.gradient(...) to evaluate grad, or return a functor which is the gradient as another function..  Lets call it math_function for now.  The behavior is things like the following:
using lambda::_1;
using lambda::_2;
math_function<double (double)> mysquare(_1 * _1); //This will generate a math_function with the finite differences as its derivative.
assert(mysquare(1) == 1);
math_function<double (double)> mysquaregrad = mysquare.gradient(); //Gets the derivative, which was generated through finite differences.  This could then be used to get the .gradient again for example.

math_function<double (double, int)> myfunc(_1 + _2);
math_function<double (double)> myboundfunc = bind(my_func, _1, 3);

//And then to create one with an analytical derivative
math_function<double (double)> mysquare(_1 * _1, 2 * _1);
assert(mysquare.gradient(1.0) == 2.0); //evaluated through the analytical derivative.
boost::function<double (double)> myboostsquare = mysquare.as_function(); //would be nice for interoperability... or something similar.

//Then we can create a class such as polynomial which subclasses math_function
polynomial_function<3> mypoly({{2, 3, 5}}); //Creates a third order polynomial.  The operator() uses fancy algorithms for polynomial evaluation.
math_function<double (double)> mypoly2 = mypoly; //Can adapt to a math_function for use in generic algorithms.