Start with replacing the indefinite integrals by inert definite integrals, so replace g2:=int(U2y,x); by g2:=Int(U2y,x=0..X); (etc) and do evalf(g2) when you have assigned values to X and k.

(1) Indefinite integration:
Your M5, M6, M7 are indefinite integrals, which cannot be evaluated numerically. For instance:

 A := Int( x^2+x, x ):  x := 0.1: evalf(A);
Error, (in evalf/int) invalid arguments

So try to do something like: M5 := Int(M4, x=0..X);
NB if you intend to evaluate all integrals numerically, use the inert form Int (not int, which tries to find an analytic primitive.

(2) Wrong definition of functions
Apparently you intend eta[1] and eta[2] to be functions of x. Then you must define these as

 eta[1]:= x -> beta1+epsilon[1]*cos(2*Pi*x);

Replace

V := select(type, [A], positive);

by

V := op( select(type, [A], positive) );

or

V := select(type, [A], positive)[];

Your equation is quite easily solved for Vb, so you get the inverse of the function you want:

VB := solve(eq,Vb);

Now you want the graph for Vb=0 to Vb=40, so you have to find the ph-values for these.
We remove the complex solutions:

ph0 := op(remove( has,[solve(VB=0)], I ));
ph40 := op(remove( has,[solve(VB=40)], I ));

Plot Vb on the horizontal axis:

plot( [VB,ph, ph=ph0..ph40], labels=["Vb","ph"] );

(this is the same as your implicitplot)
Now the derivative in, say, ph=10:

ph10 := op(remove( has,[solve(VB=10)], I )):
eval( diff( VB,ph ), ph=ph10 ): `dph/dVb(10)`:= 1/%;

f := x -> 3*x^3 - 7*x^2 + 2: #define f with arrow notation
x0 := solve( D(f)(x)=0, x );
                                14
                             0, --
                                9 
f(x0[1]); # first solution
                               2
f(x0[2]); # second solution
                              -886
                              ----
                              243 

The most straightforward way is perhaps

fsolve( 1/XX(t) = 0, t );

It will return NULL if there are no singularities or the (sequence of) t-value(s) if there are.

Did you simplify your answer?

restart;
with(LinearAlgebra):
k,i:=<sqrt(2)-1,-1,1>,<1,sqrt(2)-1,0>:
j := CrossProduct(k,i):
expand(DotProduct(k,j)),DotProduct(i,j);
                              0, 0

Use the Map procedure (capital M) from the LinearAlgebra package. Beware that the mapping is inplace, so you have to use a copy of the matrix:

restart;
A := Matrix(2,2,(i,j)->randpoly([x,y]));
B := copy(A):
LinearAlgebra:-Map( mtaylor, B, [x,y], 3 ) ;
A,B;

n := 4:
B := 'B':
LL := [[0$(n+1),1],[seq(B[i,i+1],i=-1..n),4], 
  [1,seq(B[i,i],i=0..n+1),1], [4,seq(B[i+1,i],i=0..n+1)], [1,0$(n+1)]]:
A := LinearAlgebra:-BandMatrix(LL):

Your errormessage is:

Error, unexpected `local` declaration in procedure body

This is caused by the improper use of the semicolon after the procedureheading.
And you missed a multiplication sign and a comma.

DMS17.mw

See the wikipedia page about the Cauchy–Lipschitz theorem. Write the system in the form y' = F(t,y). You can calculate the jacobian of F for the initial values which can provide you with a sufficient condition for local Lipschitz continuity.

With each iteration the number of the solutions doubles: f⁽ⁿ⁾(x)=x has 2ⁿ solutions. So for n=4 there are 16 solutions of which 12 are the zeros of a 12th degree polynomial. There is no chance that you can get explicit formulas for these roots, only numerical approximations.

restart;
f:= x -> s*x*(1-x);
h := f@@4:
S := solve( h(x)=x, {x} );

To get the values of the last 12 solutions for a certain value of s:

allvalues( eval(S[5],s=5.2) );

or, multiplied by s:

allvalues( eval( s*rhs(op(S[5])),s=5.2) );

See ?rsolve

restart;
rsolve( Y(k+1)=-10*Y(k), Y(k) ); #capital Y
                                    k
                          Y(0) (-10) 

y := unapply(%,k): #lowercase y

y(3);
                           -1000 Y(0)

Replace the curly brackets { } in the definition of h by parentheses ( ):

h := (x, y, v, w, c) -> w*(v^2+(1/3)*w*v-v-3*x-2*w-3*y)/(3-w);

restart; with(plots):
a1 := polarplot(1, t = 0 .. 6*Pi*(1/3), thickness = 2, color = black):
a4 := polarplot(1, t = 0 .. 6*Pi*(1/3),  filled = true, color = white);
a2 := polarplot(2*cos(t), t = 0 .. 6*Pi*(1/3), thickness = 2);
a3 := polarplot(2*cos(t), t = -(1/3)*Pi .. (1/3)*Pi, filled = true, thickness = 2, color=red);
display([a4,a1, a2, a3]);