@brian abraham It seems that I also still have problems. No clue why it worked several times yesterday.

I tried in Classic Worksheet Maple 14 a couple of times. There I saw no problem (using the original code). So it may be a Java problem.

@hirnyk If you ask for values (when x(0)=1, D(x)(0)=0) you get

> sol(.1);
Error, (in sol) cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

@hirnyk If you ask for values (when x(0)=1, D(x)(0)=0) you get

> sol(.1);
Error, (in sol) cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

I didn't experience any problem when trying your code.

Maple 14.01 Standard interface, Windows XP.

Algebraic:-RecursiveDensePolynomials seems to be undocumented. No help available?

Algebraic:-RecursiveDensePolynomials seems to be undocumented. No help available?

@PatrickT My reason for the type check X::name was that the procedure is not defined to handle vector input as argument number 2. Thus you should be told (in the form of an error message) if you try to pass any other object than a name.

The reason for the line

S:=select(has,indets(F(X),indexed),X);

was to allow functions containing indexed names other than the variables, e.g.

H:=X->Vector([a[8]*f(X[1],X[2]),g(X[1],X[2]),h(X[1],X[2])]);

where a[8] is just a constant.

@PatrickT My reason for the type check X::name was that the procedure is not defined to handle vector input as argument number 2. Thus you should be told (in the form of an error message) if you try to pass any other object than a name.

The reason for the line

S:=select(has,indets(F(X),indexed),X);

was to allow functions containing indexed names other than the variables, e.g.

H:=X->Vector([a[8]*f(X[1],X[2]),g(X[1],X[2]),h(X[1],X[2])]);

where a[8] is just a constant.

@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);

@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;