If I'm reading your question correctly: I don't believe there is an algorithm that, given the algorithm for your function to be integrated and appropriate initial conditions, will give an algorithm that corresponds to the integral of your original function.
However: you might wish to look into the
Chebfun project by Trefethen, Battles, Driscoll, and others. What this system does is to internally represent a function given to it as a piecewise polynomial of possibly high degree, interpolated at appropriately shifted and scaled "Chebyshev points" (roots of the Chebyshev polynomial of the first kind). The resulting
chebfun() object is then easily differentiated, integrated, or whatever other operation you might wish to do to the function. See the user guide for more details on this approach.