Question: How do I extract the data of multiple plots into a single .dat file?

Dear all,

I have 3 ODE's that I solved with the use of a for loop that cycles through different IC's. The purpose of this is so that I can compare the solutions with different IC's. I have three dependent variables (h1(z), h2(z), h3(z)), and one independent variable (z). The loop that I created looks like this:

_Omega_r := 10^(-5);
n := 10;
for i from 0 to n do:
    h10[i] := 1-0.16666667e-1*i;
    h20[i] := 1-(.7*0.16666667e-1)*i;
    h30[i]:= 3-h10[i]-h20[i];
    ans := dsolve([eq2, eq3, eq4, h1(0) = h10[i], h2(0) = h20[i], h3(0) = h30[i]], {h1(z), h2(z), h3(z)}, numeric, output = listprocedure);
    h1_sol[i] := rhs(ans[2]);
    h2_sol[i] := rhs(ans[3]);
    h3_sol[i] := rhs(ans[4]);
    print(i);
end do

Here, h10[i], h20[i] and h30[i] are my initial conditions. I use h1_sol[i], h2_sol[i] and h3_sol[i] to make plots of the solutions:

I make two different types of plots - one type has only the different solutions for h1(z) with all initial conditions:

h1plot := plot([seq(h1_sol[i](z), i = 0 .. 10)], legend = [h1(0) = h10[0], h1(0) = h10[1], h1(0) = h10[2], h1(0) = h10[3], h1(0) = h10[4], h1(0) = h10[5], h1(0) = h10[6], h1(0) = h10[7], h1(0) = h10[8], h1(0) = h10[9], h1(0) = h10[10]], title = "z versus h1(z) for different initial conditions", labels = ["z", "h1(z)"]);

And the other uses one set of initial conditions, but all three h1(z), h2(z), h3(z):

allplot1 := plot([seq(h1_sol[i](z), i = 2), seq(h2_sol[i](z), i = 2), seq(h3_sol[i](z), i = 2)], title = [typeset(h(z), " versus ", z, " (case 2)")], labels = [z, h(z)], legend = [typeset(h[1](z), " with ", h1__0 = h10[2]), typeset(h[2](z), " with ", h2__0 = h20[2]), typeset(h[3](z), " with ", h3__0 = h30[2])]);

Now, all of this works wonderfully. I managed to do everything I wanted. Except that I have been struggling figuring out how to extract the data from the plots into a .dat file into different columns. With the code that I did write, if I use it on the first plot, I end up getting 2 columns, one for z and one for h1(z) (for example), but there should be 11 columns (one for z and 10 for h1(z) with different initial conditions). How would I go about to extract the data of multiple plots into one .dat (or .xls, or .csv or whatever) file? What I did:

write := proc (plotname, filename)
    local fd; fd := fopen(filename, WRITE, TEXT);
    writedata(fd, convert(op([1, 1], plotname), matrix));
    fclose(fd);
end;

write(h1plot, `zvsh1GR.dat`)

Alternatively, if there is a way to use the solution in Maple directly without extracting first (I need to use the solutions for h1(z), h2(z), h3(z) to plot different equations that are dependent on these variables), that would be great. For example, I would like to use these three variables in an equation that uses both the derivatives of all three the variables and the variables themselves as is. For example:

And then I want to plot q(z) vs z.

I'm sorry if this is a lengthy question (or one with a simple solution) - I am still learning Maple.

My Maple file:
Solving_GR_f(R)_ne_2.maple

Maple worksheet: Solving_GR_f(R)_ne_2.mw