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

A cylinder can be regarded as the result of parallel translation of a circle:

plots[animate](plot3d,[[r*cos(phi), r*sin(phi), L], phi = 0 .. 2*Pi, r = 0 .. 3, filled = true, axes = normal, style = surface, view = [-4.7 .. 4.7, -4.7 .. 4.7, -2.7 .. 13.7], lightmodel = light4, numpoints=10000], L=0..12, frames=50);

plot3d([r*cos(phi), r*sin(phi), 12], phi = 0 .. 2*Pi, r = 0 .. 3, filled = true, axes = normal, style = surface, view = [-4.7 .. 4.7, -4.7 .. 4.7, -2.7 .. 13.7], lightmodel = light4);

Your integral can be easy calculated  numerically, if you set values  ​​n  and  h .

Example:

restart;

eta:=1000: B:=2.5: n:=1: h:=10:

evalf(Int(B*eta^(-B)*t^(B-1)*exp(-(t/eta)^B)*(t-n*h), t = n .. (n+1)*h));

                                          0.0002427195843

The global minimum is obvious without any calculation  f(0,0)=0. Initial point for finding the global maximum is easy to find from the graph:

plot3d(x^2 + y^2 + 25*(sin(x)^2+sin(y)^2), x=-2*Pi .. 2*Pi , y= -2*Pi .. 2*Pi, numpoints=3000);

Optimization[Maximize](x^2 + y^2 + 25*(sin(x)^2+sin(y)^2), initialpoint = {x =5, y=5});

                     [96.2898370467124068, [x = 4.91441834836532454, y = 4.91441834836532454]]

 Addition. The problem can be solved without any plots with  DirectSearch package:

DirectSearch[GlobalOptima](x^2 + y^2 + 25*(sin(x)^2+sin(y)^2), {x>=-2*Pi, x<=2*Pi , y>=-2*Pi,y<=2*Pi},maximize);

                           [96.2898370467124, [x = 4.9144183478452, y = -4.9144183485987], 347]

This package is available for free download link  http://www.maplesoft.com/applications/view.aspx?SID=101333