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

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