There is no attached file.

In your worksheet, you use both N and n. Did you intend for those to be the same? Does U have columns 0..m? or 1..m? or something else?

@pallav Sorry that I did not anticipate this already. Try this:

restart:
a := 2; b := .29; d := 1.85;
k:= 0;
for h from .5 by .1 to 1 do
     eq1 := x*(-b*x^2-x+1);
     eq2 := y*((a*x*x)/(b*y^2)-d-h*y);
     S := solve({eq1, eq2}, {x, y});
     SS := solve(subs(S[3], {omega^4+(h*y+x)*omega^2+h^2*x-y}), {omega});
     k:= k+1;
     tau[k] := simplify(subs(S[3], subs(SS[3], (b^2*h*y+a*x)/omega)));
end do:
convert(tau, list)[];

@pallav Sorry that I did not anticipate this already. Try this:

restart:
a := 2; b := .29; d := 1.85;
k:= 0;
for h from .5 by .1 to 1 do
     eq1 := x*(-b*x^2-x+1);
     eq2 := y*((a*x*x)/(b*y^2)-d-h*y);
     S := solve({eq1, eq2}, {x, y});
     SS := solve(subs(S[3], {omega^4+(h*y+x)*omega^2+h^2*x-y}), {omega});
     k:= k+1;
     tau[k] := simplify(subs(S[3], subs(SS[3], (b^2*h*y+a*x)/omega)));
end do:
convert(tau, list)[];

@pallav Then you need to save the taus while in the loop, and use a statement after the end of the loop to display the values. Like this:

restart:
a := 2; b := .29; d := 1.85;
for h from .5 by .1 to 1 do
     eq1 := x*(-b*x^2-x+1);
     eq2 := y*((a*x*x)/(b*y^2)-d-h*y);
     S := solve({eq1, eq2}, {x, y});
     SS := solve(subs(S[3], {omega^4+(h*y+x)*omega^2+h^2*x-y}), {omega});
     tau[h] := simplify(subs(S[3], subs(SS[3], (b^2*h*y+a*x)/omega)));
end do:
convert(tau, list)[];

@pallav Then you need to save the taus while in the loop, and use a statement after the end of the loop to display the values. Like this:

restart:
a := 2; b := .29; d := 1.85;
for h from .5 by .1 to 1 do
     eq1 := x*(-b*x^2-x+1);
     eq2 := y*((a*x*x)/(b*y^2)-d-h*y);
     S := solve({eq1, eq2}, {x, y});
     SS := solve(subs(S[3], {omega^4+(h*y+x)*omega^2+h^2*x-y}), {omega});
     tau[h] := simplify(subs(S[3], subs(SS[3], (b^2*h*y+a*x)/omega)));
end do:
convert(tau, list)[];

Could you give an example? It is not clear to me why you presented the lists three times.

If they are lists, you should enclose them in []: [a,a,a,b,c]. Otherwise, they are called sequences, which are more difficult to work with than true lists.

Do you want the rules to apply to sublists also? That is, should [b,a,a,a,b,c,e] become [b,b,b,b,c,e,e]?

@Markiyan Hirnyk That is not an example, and you know it! I ask for a solution set returned by solve which contains a trivial equation. And that is what the Asker is referring to in his followup to your Answer.

@Markiyan Hirnyk That is not an example, and you know it! I ask for a solution set returned by solve which contains a trivial equation. And that is what the Asker is referring to in his followup to your Answer.

@Markiyan Hirnyk Please kindly show me an example of a solution set returned by solve, or a similar command, which contains a trivial equation and yet represents only a finite number of solutions.

@Markiyan Hirnyk Please kindly show me an example of a solution set returned by solve, or a similar command, which contains a trivial equation and yet represents only a finite number of solutions.

Here's how you can quickly check if any solution set (or list) S contains a trivial equation of the form z = z. Do

evalb(`*`(map(lhs-rhs, S)[]) = 0);

The result will be true if and only if S contains a trivial equation.

Here's how you can quickly check if any solution set (or list) S contains a trivial equation of the form z = z. Do

evalb(`*`(map(lhs-rhs, S)[]) = 0);

The result will be true if and only if S contains a trivial equation.

@fixedeggs But dsolve itself is evaluating the derivatives by using the right sides of dsys8. There is no "built-in interpolant" for this particular step. 

But if you really want to do it your way, you just need to add one equation and one initial condition to the system for each derivative that you want to evaluate. For example,

Sys:= {diff(y(t),t) = y(t)}: IC:= {y(0)=1}:
Sys1:= {diff(y(t), t) = w(t)}: IC1:= {w(0)=1}:
Sol:= dsolve(`union`(Sys,Sys1,IC,IC1), {y(t), w(t)}, numeric);
          proc(x_rkf45_dae)  ...  end;
Sol(1);
          [t = 1., w(t) = 2.7182818283577044,  y(t) = 2.7182818283577044]

@fixedeggs But dsolve itself is evaluating the derivatives by using the right sides of dsys8. There is no "built-in interpolant" for this particular step. 

But if you really want to do it your way, you just need to add one equation and one initial condition to the system for each derivative that you want to evaluate. For example,

Sys:= {diff(y(t),t) = y(t)}: IC:= {y(0)=1}:
Sys1:= {diff(y(t), t) = w(t)}: IC1:= {w(0)=1}:
Sol:= dsolve(`union`(Sys,Sys1,IC,IC1), {y(t), w(t)}, numeric);
          proc(x_rkf45_dae)  ...  end;
Sol(1);
          [t = 1., w(t) = 2.7182818283577044,  y(t) = 2.7182818283577044]