Herd size that satisfies all the constraints, very large, and Maple can count it only approximately. Here are calculations using your link above:

w := 300426607914281713365*sqrt(609)+84129507677858393258*sqrt(7766):

n := evalf((w^4658-1/w^4658)^2/(4657*79072));

L := [7460514, 7358060, 10366482, 4149387, 4893246, 3515820, 7206360, 5439213]:

N := `+`(op(L));

N*n;

Brian, you can solve your problem by  solve  command, if  d  will be equal the least common multiple of all denominators.

restart;

unprotect(D):

L:=[B,D,W,Y,b,d,w,y]:

solve({W=5/6*B+Y, B=9/20*D+Y, D=13/42*W+Y, w=7/12*(B+b),

  b=9/20*(D+d), d=11/30*(Y+y), y=13/42*(W+w)});

assign(%):

d := ilcm(op(map(denom, L)));

L;

You can simply write:

all_plots:=[plot(sin(10*x+0.2),x=0..1, thickness=10, color=blue), plot(1-sin(10*x),x=0..1, thickness=10, color=red), plot(sin(10*x),x=0..1, thickness=10, color=green)]: 

plots:-display( all_plots );

 Done in M 12 classic.

Idea of Markiyan in above link was used.

f := piecewise(0 <= x and x <= 2, x^2, 2<x and x<3, undefined, 3 <= x and x <= 4, x^2-3, 4<x and x<5, undefined, 5 <= x and x <= 6, x^2-10):

plot(f, x = 0 .. 6, tickmarks = [default, [seq(i = i, i = 1 .. 4), 5 = "...", seq(j = j+3, j = 6 .. 13), 14 = "...", seq(k = k+10, k = 15 .. 26)]], thickness=2, discont = true, title = "Plot of  x^2");

dsolve(-(diff(lambda(s), s))-2*(diff(lambda(s), s, s))-(diff(lambda(s), s, s, s)) = 0):

subs({_C1=C[1], _C2=C[2], _C3=C[3]}, %);

The first result is also incorrect. Should be:  undefined if  -100<n<-10,  -ln|n+10|+ln|n+100| otherwise, if we consider  ln(0)=-infinity

Unfortunately, Maple weakly solves examples with parameters. Here are another two examples where the answer is incorrect:

solve(a*x = b, x);

solve(sqrt(x-a) = x, x);

Mathematica correctly solves both of these examples:

Vectors can be constructed by  plots[arrow]  command.

Example:

A := Matrix(3, 3, [1, 2, 3, 5, 1, 6, 6, 3, 9]);

V1, V2, V3:=A[..,1], A[..,2], A[..,3];

P:=plots[arrow]({V1, V2, V3}, color=red, width=[0.1, relative=false],scaling=constrained, axes=normal, orientation=[45,75]):

T:=plots[textplot3d]([[1, 5, 6, "V1"], [2, 1, 3, "V2"], [3, 6, 9, "V3"]], color=black):

plots[display](P, T, view=-1..10);

Somewhat easier to just use plot command:

f:=x->x^3+3*x^2+1:

 plot([seq([x, f(x)], x=-10..10)]);

Your procedure works correctly.

Example:

Digits := 20:

hw1(g, 0., 9):

seq(x[i], i = 0 .. 10);

 Explanation:  by default, the output of a procedure call is the last calculation.

Your spheres can be plotted by the code in single line without loops. For clarity, I reduced the number of spheres to 5:

plots[implicitplot3d]([seq(x^2+y^2+z^2 = (2*i)^2, i = 1 .. 5)], x = -10 .. 10, y = -10 .. 0, z = -10 .. 10, style = surface, numpoints = 10000, view = [-11 .. 11, -11 .. 7, -11 .. 11], scaling = constrained, axes = normal, orientation = [50, 75]);

Visualization of the surface  z=y^2+y^3+(y^3-x^2-3*x*y)^(1/4) - 5*x*y . Surface itself and its edge are built by individual  spacecurves , as  plot3d  is not completing the construction at the edges, even at high values of  numpoints . Points of global minimum are red:

F := y^2+y^3+(y^3-x^2-3*x*y)^(1/4)-5*x*y:

Sol := solve(y^3-x^2-3*x*y, x):

A := plots[spacecurve]([Sol[1], y, eval(F, x = Sol[1])], y = -3 .. 6, thickness = 2, color = blue, numpoints = 10000):

B := plots[spacecurve]([Sol[2], y, eval(F, x = Sol[2])], y = -3 .. 6, thickness = 2, color = blue, numpoints = 10000):

C1 := seq(plots[spacecurve]([x, y, F], x = Sol[1] .. Sol[2], color = blue, numpoints = 10000), y = -2.2 .. 0, .2):

C2 := seq(plots[spacecurve]([x, y, F], x = Sol[2] .. Sol[1], color = blue, numpoints = 10000), y = 0 .. 6, .2):

P1 := plottools[point]([0, 0, 0], .2, symbolsize = 12, symbol = solidsphere, color = red):

P2 := plottools[point]([4, 4, 0], .2, symbolsize = 12, symbol = solidsphere, color = red):

plots[display](P1, P2, A, B, C1, C2, view = [-10 .. 6.5, -3 .. 5.5, -4 .. 40], axes = normal, orientation = [20, 80]);

Unfortunately, the solution with Maple is rather cumbersome.

1) Consider the function  F:=y^2+y^3+(y^3-x^2-3*x*y)^(1/4) - 5*x*y :
    Construct the domain of the function  F  (green color) and find the equations of the individual branches, bounding it:

P:=plots[implicitplot](-y^3+x^2+3*x*y, x=-5..5,y=-3..5, filledregions, coloring=[green, pink],thickness=4, numpoints=10000, scaling=constrained):

T:=plots[textplot]([[-4.5, 1, "A"], [-0.2, -0.1, "B"], [1, -1.9, "E"], [3.3, -2.4, "F"], [2, -0.5, "G"], [4, 3.7, "C"]], font=[TIMES,ROMAN,18]):

plots[display](P,T);

Sol:=solve(y^3-x^2-3*x*y, x);

AB:=expand(eval(F, x=Sol[2]));  #  Values ​​of the function on individual branches

BC:=expand(eval(F, x=Sol[1]));

BEF:=expand(eval(F, x=Sol[2]));

BGE:=expand(eval(F, x=Sol[1]));

2) Obviously, the function  F  is bounded from below and at some point of its domain (or on its boundary) takes its smallest value. First, we find the critical points of  F :

Dx:=diff(F, x);

Dy:=diff(F, y);

solve({Dx=0,Dy=0});  #  The unique critical point

evalf(eval(F, %));

Since  F(0,0) = 0, then the function  F  takes its minimum value  not in found critical point, and at a point on the boundary.

3) Next, look for the minimum of the function  F  on the boundary:

minimize(AB, y=0..infinity, location);

op(simplify([minimize(BC, y=0..infinity, location)]));

minimize(BEF, y=-9/4..0, location);

minimize(BGE, y=-9/4..0, location);

Thus the final answer to the inequality  y^2+y^3+(y^3-x^2-3*x*y)^(1/4) - 5*x*y<=0 :

(x, y)=(0, 0)  and  (x, y)=(4, 4) 

You can use  subs  command in this case:

restart;

 f:=x*y+x^2+y^3:

 subs({x=x1, y=x2}, f);

                     x1*x2+x1^2+x2^3

Mathematica easily solves this inequality:

Thus we have only two solutions  (0, 0)  and  (4, 4) .

Sol1:=dsolve({diff(y(x),x)=x^2-y(x)^2, y(0)=1}, numeric):

Sol2:=dsolve({diff(z(x),x)=x^2-2*z(x)^2, z(0)=2}, numeric):

P1:=plots[odeplot](Sol1, [x,y(x)], x=0..2, thickness=2, color=red):

P2:=plots[odeplot](Sol2, [x,z(x)], x=0..2, thickness=2, color=blue):

plots[display](P1, P2);