The Taylor series methods determines the solution of an ODE by its Taylor coefficients. The theory is extremely simple but the difficulty is to determine these coefficients numerically and in an unique way. To do so automatic differentiation is used where the formulas defining the ODE are given as expression trees and are then sequentially evaluated to obtain the n-the derivative or the n-th Taylor coefficient.
I find this problem very interesting. Modern C++ methods like expression templates and operator overloading will be used and the resulting code is easy to use. Here is an example how the library should be once used:
taylor.do_step( make_vector( sigma * ( arg2 - arg1 ) , R * arg1 - arg2 - arg1 * arg3 , arg1 * arg2 - b * arg3 ) , x , t , dt );
arg1, arg2 and arg3 are placeholders for the variables x,y,z of the Lorenz system. Every entry of make_vector is an expression tree. The nodes are the operations +,-,*,/ and the leaves are arg1,arg2, arg3 and doubles. The tree is the evaluated n- times and in each evaluation the n-th derivative is builded based on the basic differentiation rules. I think this is really elegant. Existing solutions have always used some preprocessing where the expression tree are build be external programs or routines.
I have implemented a preliminary version of such a method. Surprisingly, the performance was better then hand-written Fortran code.