@Alejandro Jakubi I notice that in your (interesting!) version y and z evaluate to the local x~, although it doesn't show up in your comment.

@Alejandro Jakubi I notice that in your (interesting!) version y and z evaluate to the local x~, although it doesn't show up in your comment.

@acer

Changing eval(L) to eval(L,2) keeps the assignment of z to y:

RemoveAllAssumptions:=proc() local A,Ls,L,Lc,indx;
    Ls:=map(convert,[anames('user')],string);
    L:=ListTools:-Enumerate([anames('user')]);
    Lc:=select(hastype,eval(L,2),`local`);
    if Lc=[] then NULL
    else
       indx:=map2(op,1,Lc);
       A:=[seq(Ls[i],i=indx)];
       unassign(op(map(parse,A)));
    end if
 end proc:


@acer

Changing eval(L) to eval(L,2) keeps the assignment of z to y:

RemoveAllAssumptions:=proc() local A,Ls,L,Lc,indx;
    Ls:=map(convert,[anames('user')],string);
    L:=ListTools:-Enumerate([anames('user')]);
    Lc:=select(hastype,eval(L,2),`local`);
    if Lc=[] then NULL
    else
       indx:=map2(op,1,Lc);
       A:=[seq(Ls[i],i=indx)];
       unassign(op(map(parse,A)));
    end if
 end proc:


What is k(x)?

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

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