@goli As I mentioned earlier, you could turn your equation into a differential equation. That I have done at the end of this comment.

If you use fsolve, then add a starting value to ensure that you get the positive branch. I have chosen H=1 (it shouldn't be too crucial, which positive value you choose, but don't make it too small).

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
Y := z->fsolve(eq(z), H=1):
plot(Y, -2.2 .. 5);
implicitplot(eq(z), z = -3 .. 10,H=-10..10,gridrefine=2);
plot(fdiff(Y, [1], z), z = -2.2 .. 6, caption = "The derivative of y using fdiff on Y");
yp := implicitdiff(eq(z), H, z);
plot(eval(yp, H = 'Y(z)'), z = -2.2 .. 6, caption = "The derivative of y using implicitdiff and Y");
plot(eval((1+z)*yp/H, H = 'Y(z)'), z = -1 .. 5);

#The differential equation approach:

ode:=diff(H(z),z)=subs(H=H(z),yp);
eval(eq(0),H=1);
res:=dsolve({ode,H(0)=1},numeric,output=listprocedure);
odeplot(res,[z,H(z)],-2.22..5);
Ynum:=subs(res,H(z));
Ynum(2.345678);

@goli As I mentioned earlier, you could turn your equation into a differential equation. That I have done at the end of this comment.

If you use fsolve, then add a starting value to ensure that you get the positive branch. I have chosen H=1 (it shouldn't be too crucial, which positive value you choose, but don't make it too small).

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
Y := z->fsolve(eq(z), H=1):
plot(Y, -2.2 .. 5);
implicitplot(eq(z), z = -3 .. 10,H=-10..10,gridrefine=2);
plot(fdiff(Y, [1], z), z = -2.2 .. 6, caption = "The derivative of y using fdiff on Y");
yp := implicitdiff(eq(z), H, z);
plot(eval(yp, H = 'Y(z)'), z = -2.2 .. 6, caption = "The derivative of y using implicitdiff and Y");
plot(eval((1+z)*yp/H, H = 'Y(z)'), z = -1 .. 5);

#The differential equation approach:

ode:=diff(H(z),z)=subs(H=H(z),yp);
eval(eq(0),H=1);
res:=dsolve({ode,H(0)=1},numeric,output=listprocedure);
odeplot(res,[z,H(z)],-2.22..5);
Ynum:=subs(res,H(z));
Ynum(2.345678);

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