@goli (1) You forgot to include ode3 in the dsolve command.

(2) The right hand side of ode3 is undefined at z = 0.

(3) You require A(0) = 0, but A(z) appears in the denominator in ode3.

You can solve for H,L,R as before. Dont include A(z). The equation for A(z) you solve exactly (i. e. not numerically).

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
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)));

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

Hp,Lp,Rp:=op(subs(p,[H(z),L(z),R(z)]));
map(expand,rhs(E));
evalindets(%,specfunc(anything,Int),x->Rp(z));
AA:=eval(%,beta=Hp(0.35));
plot(AA,z=0..10);

@PatrickT

The following version has the advantage that it also accepts matrices of arbitrary dimensions mxn.

TaylorExpand := proc(F,X::name,n::posint) local d,k,X1;
d:=[op(1,F(X))][1];
X1:=op(0,F(X))(op(1,F(X)),subsop~(2=[seq(X[k],k=1..d)],op(2,F(X))),op(3..-1,F(X)));
mtaylor~(F(X),X1,n)
end proc:
TaylorExpand(M,X,2);
M2 := X-> Matrix(2,3,(i,j)->f[i,j](X[1],X[2]));
M2(X);
TaylorExpand(M2,X,2);
TaylorExpand(G,X,4);

And allowing for the number of variables to be different from the row dimension:

TaylorExpand := proc(F,X::name,n::posint) local S,nv,d,k,X1;
S:=select(has,indets(F(X),indexed),X);
nv:=max(op(op~(S)));
d:=[op(1,F(X))][1];
X1:=op(0,F(X))(op(1,F(X)),subsop~(2=[seq(X[k],k=1..nv)],op(2,F(X))),op(3..-1,F(X)));
mtaylor~(F(X),X1,n)
end proc:

H:=X->Vector([f(X[1],X[2]),g(X[1],X[2]),h(X[1],X[2])]);
TaylorExpand(H,X,4);

@PatrickT

The following version has the advantage that it also accepts matrices of arbitrary dimensions mxn.

TaylorExpand := proc(F,X::name,n::posint) local d,k,X1;
d:=[op(1,F(X))][1];
X1:=op(0,F(X))(op(1,F(X)),subsop~(2=[seq(X[k],k=1..d)],op(2,F(X))),op(3..-1,F(X)));
mtaylor~(F(X),X1,n)
end proc:
TaylorExpand(M,X,2);
M2 := X-> Matrix(2,3,(i,j)->f[i,j](X[1],X[2]));
M2(X);
TaylorExpand(M2,X,2);
TaylorExpand(G,X,4);

And allowing for the number of variables to be different from the row dimension:

TaylorExpand := proc(F,X::name,n::posint) local S,nv,d,k,X1;
S:=select(has,indets(F(X),indexed),X);
nv:=max(op(op~(S)));
d:=[op(1,F(X))][1];
X1:=op(0,F(X))(op(1,F(X)),subsop~(2=[seq(X[k],k=1..nv)],op(2,F(X))),op(3..-1,F(X)));
mtaylor~(F(X),X1,n)
end proc:

H:=X->Vector([f(X[1],X[2]),g(X[1],X[2]),h(X[1],X[2])]);
TaylorExpand(H,X,4);

@PatrickT You could do like this, where the previous definition of d has been replaced by d:=[op(1,F(X))][1];

A few other changes have been necessary.

TaylorExpand := proc(F,X::name,n::posint) local d,q,i,j,X1,X2;
d:=[op(1,F(X))][1];
if type(F(X),Matrix) and not op(1,F(X))[2]=1 then error "Matrix must have one column only" end if;
X1:=Vector(d,i->[seq(X[k],k=1..d)]);
X2:=Matrix(d,1,(i,j)->[seq(X[k],k=1..d)]);
if not type(F(X),Matrix) then mtaylor~(F(X),X1,n) else mtaylor~(F(X),X2,n) end if
end proc:
M := X-> Matrix(2,1,(i,j)->f[i,j](X[1],X[2]));

TaylorExpand(M,X,3);

I understand that an assumption is that d is also the number of variables?

@PatrickT You could do like this, where the previous definition of d has been replaced by d:=[op(1,F(X))][1];

A few other changes have been necessary.

TaylorExpand := proc(F,X::name,n::posint) local d,q,i,j,X1,X2;
d:=[op(1,F(X))][1];
if type(F(X),Matrix) and not op(1,F(X))[2]=1 then error "Matrix must have one column only" end if;
X1:=Vector(d,i->[seq(X[k],k=1..d)]);
X2:=Matrix(d,1,(i,j)->[seq(X[k],k=1..d)]);
if not type(F(X),Matrix) then mtaylor~(F(X),X1,n) else mtaylor~(F(X),X2,n) end if
end proc:
M := X-> Matrix(2,1,(i,j)->f[i,j](X[1],X[2]));

TaylorExpand(M,X,3);

I understand that an assumption is that d is also the number of variables?

@PatrickT Using one of the ways of finding the dimension of a vector given by acer you could do like this

TaylorExpand := proc(F,X::name,n::posint) local d,i;
 d:=op(1,F(X));
 mtaylor~(F(X),[[seq(X[i],i=1..d)]$d],n)
end proc;

@PatrickT Using one of the ways of finding the dimension of a vector given by acer you could do like this

TaylorExpand := proc(F,X::name,n::posint) local d,i;
 d:=op(1,F(X));
 mtaylor~(F(X),[[seq(X[i],i=1..d)]$d],n)
end proc;

As much as I hate 2d-input I appreciate the effort to make input look like mathematical expressions the way they are usually written.

Thus I appreciate very much that you can define a function f, so that f(2) returns the value of f at 2. In your comment are you somehow trying to say that the use of procedures is ineffective or slow and that the use of eval is faster or more effective?

(I should note that I too have seen students define procedures in situations where they never make use of any other input than, say x.)

As much as I hate 2d-input I appreciate the effort to make input look like mathematical expressions the way they are usually written.

Thus I appreciate very much that you can define a function f, so that f(2) returns the value of f at 2. In your comment are you somehow trying to say that the use of procedures is ineffective or slow and that the use of eval is faster or more effective?

(I should note that I too have seen students define procedures in situations where they never make use of any other input than, say x.)

Your formulation suggests that this is a homework problem or similar. Am I mistaken?

@Wang Gaoteng If a tilde ~ follows a procedure (or a function) then that procedure works elementwise on a 'container' object like a vector, an array, a list, a set, or a table.

Thus e.g.   f~([1,2,3,4]); results in [f(1), f(2), f(3), f(4)].

You could have used map of course.

map is more general, but the tilde-elementwise operation is very convenient, when container objects are involved.

@Wang Gaoteng If a tilde ~ follows a procedure (or a function) then that procedure works elementwise on a 'container' object like a vector, an array, a list, a set, or a table.

Thus e.g.   f~([1,2,3,4]); results in [f(1), f(2), f(3), f(4)].

You could have used map of course.

map is more general, but the tilde-elementwise operation is very convenient, when container objects are involved.

@Kamel Algeria

OK, what I mean is that solve doesn't make any real use of the  assumptions. `assuming` is sent to work, but solve does not in general make any use of the assumptions.

simplify(x) assuming diff(f(x),x)>0, diff(diff(f(x),x),x)<0;
`assuming`([x],[diff(f(x),x)>0, diff(diff(f(x),x),x)<0]);

@Kamel Algeria

OK, what I mean is that solve doesn't make any real use of the  assumptions. `assuming` is sent to work, but solve does not in general make any use of the assumptions.

simplify(x) assuming diff(f(x),x)>0, diff(diff(f(x),x),x)<0;
`assuming`([x],[diff(f(x),x)>0, diff(diff(f(x),x),x)<0]);

@Kamel Algeria Try this

f:=x->sin(x^3+4):
res:=solve(f(x)=0,x) assuming diff(f(x),x)>0, diff(diff(f(x),x),x)<0:
seq(eval(diff(f(x),x)>0,x=res[i]),i=1..3):
evalc([%]);

the output is

[0 < 6*2^(1/3), 0 < -3*2^(1/3)-(3*I)*2^(1/3)*sqrt(3), 0 < -3*2^(1/3)+(3*I)*2^(1/3)*sqrt(3)]

Only the first result satisfies f´(x) >0. The other two don't satisfy that requirement (and the requirement doesn't make any sense for imaginary results). So assumptions were not being used.