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

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

@Kamel Algeria In your example the result is the same if you don't use assumptions.

In the help page for solve it says (emphasis added):

"If the output of the solve command is a piecewise-defined expression, then the assuming command can be used to isolate the desired solution(s). If the output is not piecewise-defined, in particular, if the output is constant, assumptions on the independent variables may be ignored. If there are parameters in the input equations, the solve command will use those assumptions in its computations."

@Kamel Algeria In your example the result is the same if you don't use assumptions.

In the help page for solve it says (emphasis added):

"If the output of the solve command is a piecewise-defined expression, then the assuming command can be used to isolate the desired solution(s). If the output is not piecewise-defined, in particular, if the output is constant, assumptions on the independent variables may be ignored. If there are parameters in the input equations, the solve command will use those assumptions in its computations."

It would help tremendously if you would provide us with the lines of code you have used, so we don't stumble around in the dark.

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