Eigenvalues is in the module (package) LinearAlgebra. So, to use it, you need to either load the package with with(LinearAlgebra) or use the module name as a prefix. I recommend the latter. So, change Eigenvectors to LinearAlgebra:-Eigenvectors.

Assign the solution returned by solve to a variable. Like this

Sol:= solve( ...the rest of your solve command);

Then do

f:= unapply(eval(f(x), Sol), x);
g:= unapply(eval(g(x), Sol), x);

I don't believe your plot. How did you make it? I get that f(2) is about 10^7.

High-degree polynomials are notoriously unstable for numeric evaluation. You should use numeric integration. (Although, as noted by nm, the regular integrator seems to give the correct value also.) The numeric integrator will not give an answer if it cannot achieve the requested accuracy. You need to use the digits option (or the epsilon option) to make the accuracy request less than your value of Digits, which is 10 by default.

evalf(Int(f(x), x= 0..2, digits= 6));

I just picked digits=6 arbitrarily. You can make it larger, as long as Digits is somewhat larger than digits.

You must spell infinity with a lowercase i.

Numeric approximation says that your integral diverges to infinity:

Integrate(BesselJ(0,x)^2, x=0..infinity);

value(%);

evalf(%);

Two issues before answering your question:

1. You need to use `*` to represent all multiplications.
2. Once your polynomial is corrected with respect to (wrt) #1, it is already collected wrt the desired subpolynomial.

To collect wrt a subpolynomial, you need to first substitute a single name for the subpolynomial:

g:= -2-k[1]*(lambda*alpha[2]*k[2]+alpha[1]*(-2-2*k[2]+k[2]*lambda^2))+k[1]*(lambda*alpha[2]*k[2]+alpha[1]*(-2-2*k[2]+k[2]*lambda^2))*lambda*(lambda*alpha[2]*k[2]+alpha[1]*(-2-2*k[2]+k[2]*lambda^2))^2:

h:= lambda*alpha[2]*k[2]+alpha[1]*(-2-2*k[2]+k[2]*lambda^2):

subs(h= E, g);

collect(%, E);

What makes you think that Maple does not have such a command? It is called unassign.

unassign('name1', 'name2', 'name3');

If you want to see how it's done

showstat(unassign);

It's pretty simple---just 12 statements.

Start your computation over again, beginning with the command

restart;

Your problem is indeed because you entered vertical:= x=0. You should use the restart to clear that. Another way to clear the value of a variable is

vertical:= 'vertical';

However, I recommend that you get in the habit of using restart to clear mistakes.

sin(4*x)^2:
expand(%);

sort(expand(subs(sin(x)^2= 1-cos(x)^2, %)));

ErrorPlot(y, xcoords= x);

It is probably the result of assigning to an indexed name and then evaluating the unindexed form. For example,

a[1]:= 2:
eval(a);

If this isn't enough information for you to diagnose your problem, then you should post your code.

I'll answer your second question first because its solution is slightly easier, although the issue is more complex.

Your code contains a definition of a function Spec thus:

Spec:= d-> expression

where expression does not explicitly contain d although it does contain variables that do contain d. Such a construction will not make Spec a function of d. To do that, you need to use unapply:

Spec:= unapply(expression, d);

If you do this, then the computation of the max and the plot will run in under a minute.

To solve your first problem, you need to take the the max from the actual plot rather than from Spec evaluated at only integer points.

P:= plot(Spec(d), d=-15..15, ...the rest of the plot options):

#Note that evalf is not needed in plot.

A:= op([1,1], P):  #The data matrix from the plot
Smax1:= max(A[..,2]):
A[..,2]:= A[..,2]/Smax1:
P:= subsop([1,1]= A, P):
print(P);

One command handles both of your questions.

P:= a_1*x*y^2*z^3+a_2*x^4*y*z^2+a_3:
coeffs(P, [x,y,z], 'M');

M;

Note that the order of the monomials and of the coefficients is not the same as in the polynomial. However, the order of the coefficients and of the monomials are guaranteed to be the same as each other.

ode2:= convert(lhs(ode)=~rhs(ode), list)[];

print~(set_)[];

                             a = 3
                             b = 5
                             c = 9
                             

For example,

plot([sin(t), t, t= -Pi..Pi]);