@Agha  

a:=3:  b:=5:  sqrt('a-b'):  %;

@chomchom  Just change the initial condition for the first equation:

schro := {diff(psi(x), x, x)-(alpha*x^4+x^2-energy)*psi(x) = 0};

ic := {psi(3) = 0, (D(psi))(3) = 1};

Ic := [{psi(3) = 0, (D(psi))(3) = -1}, % $ 2];

E:=[1.06538, 5.74795, 8.35258]:

schro1 := [seq(subs(energy = e, alpha = 0.1, schro), e=E)];

soln1 := [seq(dsolve(schro1[i] union Ic[i], {psi(x)}, type = numeric), i=1..nops(E))];

with(plots):

display(seq(odeplot(soln1[i], [x, psi(x)], -3 .. 3, color=[red,blue,green][i]), i=1..nops(E)));

@chomchom  You can use  seq  command:

schro := {diff(psi(x), x, x)-(alpha*x^4+x^2-energy)*psi(x) = 0};

ic := {psi(3) = 0, (D(psi))(3) = 1};

E:=[1.06538, 5.74795, 8.35258]:

schro1 := [seq(subs(energy = e, alpha = 0.1, schro), e=E)];

soln1 := [seq(dsolve(schro1[i] union ic, {psi(x)}, type = numeric), i=1..nops(E))];

with(plots):

display(seq(odeplot(soln1[i], [x, psi(x)], -3 .. 3, color=[red,blue,green][i]), i=1..nops(E)));

of the initial problem everyone can see in this post .

@Carl Love  I showed only the initial step - how to integrate the equation. I think that OP will guess what to do next. If not, here's the complete solution: 

 Addition:  simplify command can be omitted,  ie just write  2*%

@brian bovril  You use combinat:-choose command two times: to select the names and separately their numerical characteristics.  In this case the order of the relevant sublists is broken. Possible workaround - combine in each sublist a name and a corresponding number:

restart;

Area:=d->Pi*(1/2*d)^2:

s := Area(15);

m := Area(20);

l := Area(25);

Obj:=(s+m+l)*(1/4);

P := [[S/2,s/2], [M/2,m/2], [L/2,l/2]];

A := combinat:-choose(P);

ans:=select(a->simplify(`+`(seq(a[i,2], i=1..nops(a)))=Obj), A);

k:=nops(ans);

pos:=seq(ListTools:-SearchAll(ans[i],A),i=1..k);

@Carl Love

If we can not disprove a certain statement, then we consider it a true:

evalb(a<>x); 

                                         true

@emendes  Test1  procedure generalizes  Test  procedure to any number of variables ( f  is the list of expressions,  v  is the list of variables):

Test1:=(f, v)->`if`(`or`(seq(depends(f[i], v[i]) and not depends(f[i], {v[]} minus {v[i]}), i=1..nops(v))) , false, true):

I did not check your options (no free time).

Carl, your version can give a result that is different from mine:

f:=[a,b,c]: v:=[x,y,z]:

Test1(f, v);

valid(f, v);

                                                true

                                                false

I think OP will decide which version is more suitable to him.

@vv  I think the reason is that Maple is plotting  the surface by the values of the function on a rectangular grid. The domain of the function may be different from this rectangular region and many points do not lie on the border. Here is a simple example, where it is clearly visible:

S:=y^(1/3);

plot3d(S, x=0..1, y=0..1, grid=[30,30], axes=normal, view=[0..1,0..1,0..1], orientation=[-25,70]);  # OK

Slightly change the range for  y  variable:

plot3d(S, x=0..1,y=-0.1..1 , grid=[30,30], axes=normal, view=[0..1,0..1,0..1], orientation=[-25,70]);  # Not very good plotting because the grid's points do not lie on the  x-axis

   Change the grid in order to its points again lie on the  x-axis:

plot3d(S,  x=0..1, y=-0.1..1, grid=[23,23], axes=normal, view=[0..1,0..1,0..1], orientation=[-25,70]);  # Now OK again

@Bendesarts  I do not understand. You wrote "I have tried but it seems that no simplifications is made with this function simplify(,size)".   Compare the lengths of expressions  condition1  and  condition2 . The latter is about 5 times shorter:

Download simplification_condition_de_compatibilité_new.mw

@diedrotn  I do not know any other way to solve your equation.

You have to give a detailed description of the problem as a whole. Maybe then someone will be able to help you.

@taro   I don't know. I never used the syntax  `convert/string`(...)  , only  convert(..., string)

@vv  Thank you.

You wrote "But what would be the simplest workaround for g(x) ?" I don't know. I think that the use of a special symbol like  _x   - is a good idea.

@vv  No, I haven't any references.

@vv  Thank you for your interest. I think that your conjectured theorem is true. Moreover, it can be reinforced by any number of rotating circles  n>=3 . Here is an example of animation with 4 circles. Here should be  x<=3-2*sqrt(2)

Circles1:=proc(x)

local OO, O1, O2, O3, O2x, O2y, OT, dist, Sys, Sol, sol, y, u, z, P;

uses plottools, plots;

OO:=[0,0]: O1:=[x+y,0]: O2:=[O2x,O2y]: O3:=[-x-z,0]: OT:=[x+2*y-1,0]:

dist:=(A,B)->sqrt((B[1]-A[1])^2+(B[2]-A[2])^2):

Sys:={dist(O1,O2)^2=(y+u)^2, dist(OO,O2)^2=(x+u)^2, dist(O2,O3)^2=(u+z)^2,   x+y+z=1, dist(O2,OT)^2=(1-u)^2};

Sol:=op~([allvalues([solve(Sys)])]);

sol:=select(i->is(eval(convert([y>0,u>0,z>0,O2y>0,x<=y,u<=y,z<=u],`and`),i)), Sol)[];

assign(sol);

O1:=[x+y,0]: O2:=[O2x,O2y]: O3:=[-x-z,0]: OT:=[x+2*y-1,0]:

P:=proc(phi)

local eq, r1, r, R, Ot, i, S, s, t, P1, P2;

uses plots,plottools;

eq:=1-dist([r*cos(s),r*sin(s)],OT)=r-x;

r1:=solve(eq,r);

r:=eval(r1,s=phi);

R[1]:=evalf(r-x);

Ot[1]:=evalf([r*cos(phi),r*sin(phi)]);

for i from 2 to 4 do

S:=[solve({1-dist(OT,[s,t])=dist(Ot[i-1],[s,t])-R[i-1], 1-dist(OT,[s,t])=dist(OO,[s,t])-x})];

P1:=eval([s,t],S[1]); P2:=eval([s,t],S[2]);

Ot[i]:=`if`(evalf(Ot[i-1][1]*P1[2]-Ot[i-1][2]*P1[1])>0,P1,P2);

R[i]:=dist(Ot[i],OO)-x;

od;

display(circle([x+2*y-1,0],1, color=blue,thickness=4), circle(OO,x, color=red,thickness=4), seq(circle(Ot[k],R[k], thickness=3),k=1..4), scaling=constrained, axes=none);

end proc:

animate(P,[phi], phi=0..Pi, frames=60);

end proc:

Circles1(0.15);