How to integrate the contents of a vector using an adaptive quadrature routine

Question

I have a function which requires the return type to be a container. The problem is that I need to integrate the contents of the container as efficiently as possible and was hoping to use adaptive Gauss-Kronrod integration or something equally efficient (or better).

I was hoping to use the GNU Scientific Library quadrature routine qags, but it returns the result with type double. As things stand I think I'm probably going to have to re-write sections of the GSL quadrature routines to convert the return type to std vector, but that would be a rather lengthy and possibly error-prone detour. I was hoping somebody might recommend a better solution!

Below is an example of the kind of function I'm trying to integrate over, presently employing the basic trapezoidal rule but indicating where I would instead prefer to implement the GSL routine gaq. whilst it is much simpler than my actual problem, it demonstrates that each element that is put in the vector is calculated from the previous result, hence the requirement of the container.

#include <iostream>
#include <vector>
#include <gsl/gsl_integration.h>
using namespace std;

vector<double> f(double E, int N) {
    vector<double> result;
    result.reserve(N);
    double x = E;
    for (int it=0; it < N; ++it){
        result.push_back(x);
        x *= x;
    }
    return result;
}

vector<double> f_qag(double E, void * params) {
    int N = *(int *) params;
    vector<double> result;
    result.reserve(N);
    double x = E;
    for (int it=0; it < N; ++it){
        result.push_back(x);
        x *= x;
    }
    return result;
}


int main(){
    vector<double> result;
    vector<double> integrate;
    int N = 20;
    result.reserve(N);
    integrate.reserve(N);
    for (int i = 0; i < N; i++)
        result[i] = 0.;
    double end = 1.0;
    double start = -1.0;

    // I would like to use qag here
    /* gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000); */
    /* vector<double> error; */
    /* gsl_function F; */
    /* F.function = &f_qag; */
    /* F.params = &N; */
    /* gsl_integration_qag (&F, start, end, 0, 1e-7, 1000, 1, w, &result, &error); */

    // instead of resorting to the trapezoidal rule here
    double E;
    const int n = 1000;
    double factor = (end - start)/(n*1.);
    for (int k=0; k<n+1; k++) {
        E = start + k*factor;
        integrate = f(E, N);
        for (int i = 0; i < N; i++){
            if ((k==0)||(k==n))
                result[i] += 0.5*factor*integrate[i];
            else 
                result[i] += factor*integrate[i];
        }
    }   

    for (int i = 0; i < N; i++)
        cout<<i<<" "<<result[i]<<endl;

    return 0;
}

I intend to check for convergence against the pseudo-result conv;

double conv = 0.;
for (int i = 0; i < N; i++)
    conv += abs(integrate[i]);