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

Examples of two variants.

With the for loop:

a[0]:=t->0: a[1]:=t->t:

N:=5:

for n from 2 to N do

a[n]:=unapply(a[n-1](t)^2-3*a[n-1](t)*a[n-2](t)+a[n-2](t), t):

od:

sort(simplify(a[N](t)));

t^16-12*t^15+48*t^14-50*t^13-111*t^12+225*t^11+17*t^10-174*t^9+129*t^8+17*t^7-81*t^6+24*t^5+14*t^4-9*t^3+t

With the recursive procedure:

a:=proc(N)

local f;

option remember;

a(0):=t->0; a(1):=t->t;

f:=t->a(N-1)(t)^2-3*a(N-1)(t)*a(N-2)(t)+a(N-2)(t);

unapply(sort(simplify(f(t))), t);

end proc:

a(5)(t);

plot(%, t=-1..1,-2..4);

restart;

f:=proc(n)

option remember;

unapply(f(0)(t)+(1/omega)*int((omega^2*f(n-1)(tau)-f(n-1)(tau)^3)*sin(omega*(t-tau)),tau=0..t), t);

end proc:

f(0):=t->a*cos(omega*t):

Expression for function  is very quickly complicating with increasing  n. For example, if n = 3, then time is about 3.5 minutes. I did not wait for the end of the calculations for n = 4.

Example for  n=2:

f(2);

Writing  s := solve({eq[1],..,eq[6]},{x1,..,x6});  is erroneous.

Possible variant

s := solve({seq(eq[i], i=1..nops(A))}, {op(b)});

restart:

A:=Matrix([[1,2,3],[4,8,6],[7,8,9]]):

X:=Vector(3,symbol=x):

F:=Vector(3,[8,8,9]):

solve({seq((A.X)[i]=F[i], i=1..3)});

C3 := <3*cos(t), 3*sin(t), s>:

P := Matrix(3, 3, [[-(2/5)*sqrt(5), -(2/15)*sqrt(5),  -1/3], [0, (1/3)*sqrt(5),-2/3], [(1/5)*sqrt(5), -(4/15)*sqrt(5), -2/3]]):

V := <-2, 10, 3>;

C4 := P.C3;

Rightcylinder1 := plot3d([C4[1]+V[1], C4[2]+V[2], C4[3]+V[3]], t = 0 .. 2*Pi, s = 0 .. 12, numpoints = 6500, thickness = 1, style = surface, color = red):

Plane1 := plots[implicitplot3d](x-V[1]+2*(y-V[2])+2*(z-V[3]) = 0, x = -20 .. 20, y = -20 .. 20, z = -20 .. 20, style = wireframe, numpoints = 10000):

plots[display](Plane1, Rightcylinder1, axes = normal, orientation = [-55, 75], view = [-10 .. 15, -15 .. 20, -10 .. 15], scaling = constrained);

First of all it is necessary to introduce a measure of closeness of two functions. The easiest way to use the quadratic metric  distance(f, g)=sqrt(int((f(t)-g(t))^2, t=0..a))

Your example (approximation  f(t)=exp(t)  for N=4,  M=3,  a=1):

HybrFunc(4, 3, 1):

g(t):=add(add(c[n,m]*b[n,m](t), m=0..2), n=1..4);

Optimization[Minimize](sqrt(int((exp(t)-g(t))^2, t=0.. 1)));

Visualization:

assign(op(%[2]));

plot([exp(t),g(t)], t=-1..2, 0..5, color=[blue,red],thickness=[1,3],discont, scaling=constrained);

If we define the values ​​of all parameters, there is no problem.

Example:

lambda[v] := 3: lambda[t] := 1: q[p] := 2: t[c] := 4: S[di] := 3: h := 10:

int(int(lambda[v]*lambda[t]*exp(-lambda[v]*v-lambda[t]*t), v = (1/2*(q[p]+q[p]*t[c]*t+2*S[di]*h*t))/(h*t) .. infinity), t = 0 .. infinity);

evalf(%);

With arbitrary parameters Maple does not evaluate even more simple integral:

int(int(lambda[v]*lambda[t]*exp(-lambda[v]*v-lambda[t]*t), v = 0 .. infinity), t = 0 .. infinity)  assuming lambda[v] > 0, lambda[t] > 0;

For greater clarity, I made to your code some minor changes:

1) In order to the cylinder be above the plane  x+2*y+2*z = 0  but not below, the direction vector  w3  is reversed.

2) Changed the color of the cylinder and the grid on the plane became more frequent.

Corrected lines of the code:

w3:=<1/3, 2/3, 2/3>:                  

Rightcylinder := plot3d([C2[1], C2[2], C2[3]], t = 0 .. 2*Pi, s = 0 .. 12, numpoints = 4000, thickness = 2, scaling = constrained, color = "LightBlue"):

Plane := implicitplot3d(x+2*y+2*z = 0, x = -20 .. 20, y = -20 .. 20, z = -20 .. 20, style = wireframe, numpoints = 10000):

a := s*w3:

Axis := spacecurve([a[1], a[2], a[3]], s = 0 .. 12, color = red, thickness = 2):

 display(Plane, Rightcylinder, Axis,  orientation=[45, 60], axes = normal, view = [-10 .. 10, -5 .. 10, -5 .. 10]);

Output is a list of lists. Each sublist consists of a number of the relevant member of the expression and the list of the member and of the coefficient in front of it:

> Ex:=(1/20)*t[1]^2*(diff(phi(x), x, x))^2*((t[1]+t[2])^5-t[1]^5)/(t[2]^2*E[T])+(1/8)*t[1]^2*(diff(phi(x), x, x))^2*(-2*t[1]-2*t[2])*((t[1]+t[2])^4-t[1]^4)/(t[2]^2*E[T])+(1/3*((1/4)*t[1]^2*(diff(phi(x), x, x))^2*(t[1]+t[2])^2/(t[2]^2*E[T])+(1/4)*t[1]^2*(diff(phi(x), x, x))^2*(-2*t[1]-2*t[2])^2/(t[2]^2*E[T])+(1/2)*t[1]*(diff(phi(x), x, x))*((-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*phi(x)/t[2]+(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*eta(x)/t[2]+(1/2)*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/(E[T]*t[2])+2*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*psi(x)/t[2])/t[2]-t[1]*(diff(psi(x), x))*(-(cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))/t[2]-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))/t[2])/t[2]+(1/2)*t[1]^2*eta(x)*(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*(diff(phi(x), x, x))/t[2]^2+(1/2)*t[1]^2*phi(x)*(-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*(diff(phi(x), x, x))/t[2]^2-t[1]*(diff(phi(x), x))*(-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))/t[2]-(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))/t[2])/t[2]+t[1]^2*psi(x)*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*(diff(phi(x), x, x))/t[2]^2))*((t[1]+t[2])^3-t[1]^3)+(1/2*(t[1]^2*psi(x)*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*(diff(phi(x), x, x))*(-2*t[1]-2*t[2])/t[2]^2+t[1]*(diff(phi(x), x))*(t[1]+t[2])*(-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))/t[2]-(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))/t[2])/t[2]-t[1]*(diff(phi(x), x))*(cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2])/t[2]+t[1]*(diff(psi(x), x))*(t[1]+t[2])*(-(cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))/t[2]-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))/t[2])/t[2]-t[1]*(diff(psi(x), x))*((cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2])/t[2]+(1/2)*t[1]^2*eta(x)*(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*(diff(phi(x), x, x))*(-2*t[1]-2*t[2])/t[2]^2+(1/2)*t[1]^2*phi(x)*(-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*(diff(phi(x), x, x))*(-2*t[1]-2*t[2])/t[2]^2+(1/4)*t[1]^2*(diff(phi(x), x, x))^2*(t[1]+t[2])^2*(-2*t[1]-2*t[2])/(t[2]^2*E[T])+(1/2)*t[1]*(diff(phi(x), x, x))*(-2*t[1]-2*t[2])*((-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*phi(x)/t[2]+(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*eta(x)/t[2]+(1/2)*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/(E[T]*t[2])+2*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*psi(x)/t[2])/t[2]))*((t[1]+t[2])^2-t[1]^2)+t[1]*phi(x)*((cos(a)^4/E[A]+cos(a)^2*sin(a)^2*(-2*nu[A]/E[A]+1/G[A])+sin(a)^4/E[T])*t[1]*phi(x)/t[2]+(sin(a)^2*cos(a)^2*(1/E[A]+1/E[T]-1/G[A])-(sin(a)^4+cos(a)^4)*nu[A]/E[A])*t[1]*eta(x)/t[2]+(1/2)*(-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/t[2]+sin(a)*cos(a)*(cos(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+sin(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*psi(x)/t[2])+t[1]*(diff(phi(x), x))*(t[1]+t[2])*(cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2])+(1/2)*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2*((-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*phi(x)/t[2]+(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*eta(x)/t[2]+(1/2)*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/(E[T]*t[2])+2*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*psi(x)/t[2])+t[1]*eta(x)*((sin(a)^2*cos(a)^2*(1/E[A]+1/E[T]-1/G[A])-(sin(a)^4+cos(a)^4)*nu[A]/E[A])*t[1]*phi(x)/t[2]+(sin(a)^4/E[A]+cos(a)^2*sin(a)^2*(-2*nu[A]/E[A]+1/G[A])+cos(a)^4/E[T])*t[1]*eta(x)/t[2]+(1/2)*(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/t[2]+sin(a)*cos(a)*(sin(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+cos(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*psi(x)/t[2])+t[1]*psi(x)*(sin(a)*cos(a)*(cos(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+sin(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*phi(x)/t[2]+sin(a)*cos(a)*(sin(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+cos(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*eta(x)/t[2]+cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/t[2]+(4*sin(a)^2*cos(a)^2*(1/E[A]+2*nu[A]/E[A]+1/E[T])+(sin(a)^2-cos(a)^2)^2/G[A])*t[1]*psi(x)/t[2])+t[1]*(diff(psi(x), x))*(t[1]+t[2])*((cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2]):

Ex1:=expand(Ex);

nops(Ex1);

L:=[seq([i,[op(i,Ex1),`if`(type(op(1,op(i,Ex1)),realcons),op(1,op(i,Ex1)), 1)]], i=1..nops(Ex1))]:

L[1..20];