Thank you but I really did it. I tried maxfun=500000 and also maxfun=1000000 and more, many times but the error does not solve.  

thanks a lot

thanks a lot

thank you. I did it myself. can you propose another way for it?

thank you. I did it myself. can you propose another way for it?

@hirnyk Thank you very much. It works very nice.

@hirnyk Thank you very much. It works very nice.

amin.mw Thanks for your reply. It's not exactly what I want. So, I've attached my file.

I can obtain the value of "L" in terms of any "t" (you see). But I need the values of "L" in terms of "zi". I've written the relation between "t" and "zi", as well.

Thanks again!

amin.mw Thanks for your reply. It's not exactly what I want. So, I've attached my file.

I can obtain the value of "L" in terms of any "t" (you see). But I need the values of "L" in terms of "zi". I've written the relation between "t" and "zi", as well.

Thanks again!

@Robert Israel Thanks Mr. Israel! You are right. I had not written an array(I was wrong), but I have an array.

@Robert Israel Thanks Mr. Israel! You are right. I had not written an array(I was wrong), but I have an array.

@DuncanA Dear my friend!

One new related question. Consider I want to identify " m "s and " n "s, related to those values of alpha (the third column) which are smaller than for example " 5 ". I need a new array shows these " m "s and " n "s. How can I do that?

thanks a lot!

@DuncanA Dear my friend!

One new related question. Consider I want to identify " m "s and " n "s, related to those values of alpha (the third column) which are smaller than for example " 5 ". I need a new array shows these " m "s and " n "s. How can I do that?

thanks a lot!

Thanks for your reply. I think you know maple very well. I want to do your notation for 4 following ODEs:

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
Y := z->if not type(z,numeric) then 'procname(z)' else fsolve(eq(z), H=1) end if:

yp := implicitdiff(eq(z), H, z);

ode:=diff(H(z),z)=subs(H=H(z),yp);

ode1:=diff(L(z),z)=L(z)/(1+z)+(1+z)/H(z);

ode2:=diff(R(z),z)=1/H(z);

ode3:=diff(A(z),z)=(-2*A(z))/(3*z)+(beta^(1/3))/(3*z*H(z)*(A(z)^(1/2)));

p := dsolve({ode,ode1,ode2,L(0)=0,H(0)=1,R(0)=0,A(0)=0}, numeric, output=listprocedure);

where "beta" is a constant. This doesn't work because singularity at "z=0" for both "z" and "A(z)", I think. So, how can I do this?

Also, if "beta = H(0.35)" how should I rewrite my ODEs?

thank you very much!  

Thanks for your reply. I think you know maple very well. I want to do your notation for 4 following ODEs:

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
Y := z->if not type(z,numeric) then 'procname(z)' else fsolve(eq(z), H=1) end if:

yp := implicitdiff(eq(z), H, z);

ode:=diff(H(z),z)=subs(H=H(z),yp);

ode1:=diff(L(z),z)=L(z)/(1+z)+(1+z)/H(z);

ode2:=diff(R(z),z)=1/H(z);

ode3:=diff(A(z),z)=(-2*A(z))/(3*z)+(beta^(1/3))/(3*z*H(z)*(A(z)^(1/2)));

p := dsolve({ode,ode1,ode2,L(0)=0,H(0)=1,R(0)=0,A(0)=0}, numeric, output=listprocedure);

where "beta" is a constant. This doesn't work because singularity at "z=0" for both "z" and "A(z)", I think. So, how can I do this?

Also, if "beta = H(0.35)" how should I rewrite my ODEs?

thank you very much!