@goli Instead of me guessing what you are doing, I suggest that you bring the whole code beginning with a restart to indicate that nothing has been defined before.

@goli Instead of me guessing what you are doing, I suggest that you bring the whole code beginning with a restart to indicate that nothing has been defined before.

I cannot reproduce your result of 858.000000000000114.

I get the integer 858 (not 858.) as I expected. I'm using Maple 14, worksheet interface, Maple input (of course).

If I change the division by the integer 1000 to division by the float 1000. I get 858.000000000000.

@goli Have you tried it yourself? Did you run into any problems? If so, what were they?

@goli Have you tried it yourself? Did you run into any problems? If so, what were they?

@goli Surely the positive solution is used also for y'(x).

Try this plot, which shows the function and its derivative.

plot(['Y(x)',eval(yp,y='Y(x)')],x=-3..1.5,caption="y (i.e. Y, red) and the derivative of y (using implicitdiff and Y, blue)",color=[red,blue]);

Notice that

factor(yp);

returns .

Thus yp has two zeros, nicely confirmed by the plot.

Another way of getting to your numerical solution is to turn your equation into a differential equation with an appropriate initial condition, which could be y(-1) = 1, since x=-1 and y=1 satisfies eq(x).

ode:=diff(y(x),x)=subs(y=y(x),yp);

res:=dsolve({ode,y(-1)=1},numeric,output=listprocedure);
Yode:=subs(res,y(x));
plot(Yode,-3..1.5);

If you only want to plot y and y', then you need only do

plots:-odeplot(res,[[x,y(x)],[x,rhs(ode)]],-3..1.5);

@goli Surely the positive solution is used also for y'(x).

Try this plot, which shows the function and its derivative.

plot(['Y(x)',eval(yp,y='Y(x)')],x=-3..1.5,caption="y (i.e. Y, red) and the derivative of y (using implicitdiff and Y, blue)",color=[red,blue]);

Notice that

factor(yp);

returns .

Thus yp has two zeros, nicely confirmed by the plot.

Another way of getting to your numerical solution is to turn your equation into a differential equation with an appropriate initial condition, which could be y(-1) = 1, since x=-1 and y=1 satisfies eq(x).

ode:=diff(y(x),x)=subs(y=y(x),yp);

res:=dsolve({ode,y(-1)=1},numeric,output=listprocedure);
Yode:=subs(res,y(x));
plot(Yode,-3..1.5);

If you only want to plot y and y', then you need only do

plots:-odeplot(res,[[x,y(x)],[x,rhs(ode)]],-3..1.5);

@goli I think you are using Robert Israels's reformulation of the equation

eq2:= x->(y^2 - (1+x)^3 - (1+x)^4)^5 - y=0;

which for y>0 is equivalent to eq. However, for general y they are not equivalent.

The following does give two curves, yes, one above the x-axis, and one below.

plots:-implicitplot(eq2(x), x=-10..10,y=-10..10,grid=[100,100],gridrefine=2);

Using eq you only get a curve above the x-axis.

I thought you were interested in the one above?

@goli I think you are using Robert Israels's reformulation of the equation

eq2:= x->(y^2 - (1+x)^3 - (1+x)^4)^5 - y=0;

which for y>0 is equivalent to eq. However, for general y they are not equivalent.

The following does give two curves, yes, one above the x-axis, and one below.

plots:-implicitplot(eq2(x), x=-10..10,y=-10..10,grid=[100,100],gridrefine=2);

Using eq you only get a curve above the x-axis.

I thought you were interested in the one above?

@goli Here are the different commands.

restart;
eq:=x->y^(-1/5)*(y^2-(1+x)^3-(1+x)^4)=1;
Y:=x->fsolve(eq(x),y);
plot(Y,-10..10);
plots:-implicitplot(eq(x), x=-10..10,y=0..120,grid=[150,150]);
plot(fdiff(Y,[1],x),x=-2..2,caption="The derivative of y using fdiff on Y");
yp:=implicitdiff(eq(x),y,x);
plot(eval(yp,y='Y(x)'),x=-2..2,caption="The derivative of y using implicitdiff and Y");

@goli Here are the different commands.

restart;
eq:=x->y^(-1/5)*(y^2-(1+x)^3-(1+x)^4)=1;
Y:=x->fsolve(eq(x),y);
plot(Y,-10..10);
plots:-implicitplot(eq(x), x=-10..10,y=0..120,grid=[150,150]);
plot(fdiff(Y,[1],x),x=-2..2,caption="The derivative of y using fdiff on Y");
yp:=implicitdiff(eq(x),y,x);
plot(eval(yp,y='Y(x)'),x=-2..2,caption="The derivative of y using implicitdiff and Y");

After using combine you can use

convert(%,phaseamp,t);

After using combine you can use

convert(%,phaseamp,t);

My Maple 14 uses HeunG, no HG appears.

@GeorgesL In any case, it is not a bug. It is a decision that was made. What do you think 0^k should evaluate to, when k is unassigned?