I found that if I remove "stiff=true" it takes less time. Thank you!

@Robert Israel I used it and it works well. Thanks a lot.

@Robert Israel I used it and it works well. Thanks a lot.

@Axel Vogt I couldn't copy and paste my sheet in math notation, so first I exchanged it to classical sheet by copy and paste and then to this window. 

Thanks for your reply. For example after solving an ODE system I want to obtain the result of an integral as follows:

> b := -95; `ε` := 1/4; h := .69; n := 9.2; om := .7:

> ode1 := diff(Omega(z), z)+(Omega(z)*(1-Omega(z))*(3-2*sqrt(Omega(z))*(1+z)/n)+Omega(z)*(1-Omega(z))*`ε`*b*(1+z)*(diff(phi(z), z)))/(1+z) = 0:

> ode2 := diff(phi(z), z)+sqrt(2*(1+z)*sqrt(Omega(z))/(3*n)+(1-Omega(z))*(1+z)*b*`ε`*(diff(phi(z), z))/(3*Omega(z)))/((1+z).H(z)) = 0:

> a := solve(ode2, diff(phi(z), z)):

> ode22 := diff(phi(z), z) = a[1]:

> ode3 := diff(H(z), z)+H(z)*((3/2)*Omega(z)-(1+z)*Omega(z)^(3/2)/n+(1/2)*(1-Omega(z))*(1+z)*b*`ε`*(diff(phi(z), z))-3/2)/(1+z) = 0:

> sys := {ode1, ode22, ode3}:

   ics := {H(0) = .69, phi(0) = 0, Omega(0) = .7}:

   sol := dsolve(`union`(sys, ics), numeric, output = listprocedure, stiff = true): 

> HH := subs(sol, H(z));

The integral I want to solve is:

> R := h*evalf(sqrt(1-om)*(int(1/HH(z), z = 0 .. 1091.3)));

Thanks for your reply. For example after solving an ODE system I want to obtain the result of an integral as follows:

> b := -95; `ε` := 1/4; h := .69; n := 9.2; om := .7:

> ode1 := diff(Omega(z), z)+(Omega(z)*(1-Omega(z))*(3-2*sqrt(Omega(z))*(1+z)/n)+Omega(z)*(1-Omega(z))*`ε`*b*(1+z)*(diff(phi(z), z)))/(1+z) = 0:

> ode2 := diff(phi(z), z)+sqrt(2*(1+z)*sqrt(Omega(z))/(3*n)+(1-Omega(z))*(1+z)*b*`ε`*(diff(phi(z), z))/(3*Omega(z)))/((1+z).H(z)) = 0:

> a := solve(ode2, diff(phi(z), z)):

> ode22 := diff(phi(z), z) = a[1]:

> ode3 := diff(H(z), z)+H(z)*((3/2)*Omega(z)-(1+z)*Omega(z)^(3/2)/n+(1/2)*(1-Omega(z))*(1+z)*b*`ε`*(diff(phi(z), z))-3/2)/(1+z) = 0:

> sys := {ode1, ode22, ode3}:

   ics := {H(0) = .69, phi(0) = 0, Omega(0) = .7}:

   sol := dsolve(`union`(sys, ics), numeric, output = listprocedure, stiff = true): 

> HH := subs(sol, H(z));

The integral I want to solve is:

> R := h*evalf(sqrt(1-om)*(int(1/HH(z), z = 0 .. 1091.3)));

Thanks PatrickT 1503. I'll try to add another ODE. I think that's a good idea.

Thanks PatrickT 1503. I'll try to add another ODE. I think that's a good idea.

Thanks, but I don't want maple shows to me all the numbers 1, 2, ... in a column. I want to print for example 2 instead of 1. What should I do now?

Thanks, but I don't want maple shows to me all the numbers 1, 2, ... in a column. I want to print for example 2 instead of 1. What should I do now?

@Preben Alsholm Dear Preben. I think there are some mistakes with maple (or maybe with me). I wanted to do all of the previous equations and integrals for a large set of initial values. So, I used loop for 14641 set of initial values as:

for m from 0.25 by 0.01 to 0.35 do       

for n from -0.35 by 0.01 to -0.25 do 

for k from -0.05 by 0.01 to 0.05 do

for h from 0.65 by 0.01 to 0.75 do

Apparently, maple did all of the things that I wanted and wrote my results data in a text file. But when I limited my ranges for example for 81 set of initial values as:

for m from 0.33 by 0.01 to 0.35 do

for n from -0.36 by 0.01 to -0.34 do 

for k from 0.03 by 0.01 to 0.05 do

for h from 0.68 by 0.01 to 0.7 do

I found errors and when I checked, I thought the problem is that for some of special values maple can not solve the equations. Then, I thought maybe maple has given me wrong data for the large set. So, I asked you and you guided me to solve the problem. But when I use your notation for my large set again, maple shows exactly the previous results. So, my question is how maple did all the equations and integrals for my large set at the first time (before asking you) though it could not solve the integral for some special cases of initial values? (I should notice that your notation works for the small set, as well) 

@Preben Alsholm Dear Preben. I think there are some mistakes with maple (or maybe with me). I wanted to do all of the previous equations and integrals for a large set of initial values. So, I used loop for 14641 set of initial values as:

for m from 0.25 by 0.01 to 0.35 do       

for n from -0.35 by 0.01 to -0.25 do 

for k from -0.05 by 0.01 to 0.05 do

for h from 0.65 by 0.01 to 0.75 do

Apparently, maple did all of the things that I wanted and wrote my results data in a text file. But when I limited my ranges for example for 81 set of initial values as:

for m from 0.33 by 0.01 to 0.35 do

for n from -0.36 by 0.01 to -0.34 do 

for k from 0.03 by 0.01 to 0.05 do

for h from 0.68 by 0.01 to 0.7 do

I found errors and when I checked, I thought the problem is that for some of special values maple can not solve the equations. Then, I thought maybe maple has given me wrong data for the large set. So, I asked you and you guided me to solve the problem. But when I use your notation for my large set again, maple shows exactly the previous results. So, my question is how maple did all the equations and integrals for my large set at the first time (before asking you) though it could not solve the integral for some special cases of initial values? (I should notice that your notation works for the small set, as well) 

@Preben Alsholm Dear Preben! Thank you very much. Your notation works well. That's very nice. good luck.

@Preben Alsholm Dear Preben! Thank you very much. Your notation works well. That's very nice. good luck.

@Preben Alsholm Thanks a lot. But I have some new problems, now. After, solving this equation we want to solve another equation (ODE), that with your new version of fsolve, maple does not solve it. The ODE is:

> l := dsolve({D(L)(z) = L(z)/(1+z)+(1+z)*h/Y(z),L(0)=0}, type=numeric, range=0..2,known=Y);

Also, my final goal to find Y for large values of z was to determne some integrals of Y(z), where I couldn't get their results as a number. This problem exists yet. For example:

> A1 := ((0.35/Y(0.35))*(int(1/Y(z), z = 0 .. 0.35))^2)^(1/3);              and:

> R := h*evalf(sqrt(m)*(int(1/Y(z), z = 0 .. 1091.3)));

For the first integral I could find result before, but with your new reply I can't. What's the problem?

Also, I have a question. Is the starting value in fsolve shows the value of H at z=0? If it is right I must set H=h, for my physical concept and then what does your last reply mean?

Thank you very much.