Based on the information that you have provided, this is all that we can do:

If that's not what you want, then you will have to be more specific about what you are looking for.

Download mw.mw

Here are two different ways to achieve that.

Method 1: Use the EulerLagrange function in the VariationalCalculus package:

VariationalCalculus:-EulerLagrange(L, t, [ x(t), y(t), alpha(t), phi1(t), phi2(t) ] );
remove(type, %, equation);

Method 2: Load the Physics package:

with(Physics):

That redefines diff so that your code should work unless you have other errors in it.

There are other ways as well but the two that I have shown here are the quickest. 

In the future, when you ask a question, include your actual Maple worksheet.  Showing just a picture, as you have done, is not very helpful since I cannot execute that code to test my suggestions before replying.

The polygonplot3d help page describes three different ways of invoking that function.  Their descriptions are somewhat difficult to decipher, so here is an expanded description.

In the first method, we pass the vertex coordinates to polygonplot3d as a single list of lists, each internal list providing the coordinates of one vertex of the polygon:

plots:-polygonplot3d( [ [x1,y1,z1], [x2,y2,z2], ..., [xn,yn,zn] ]);

In the second method we pass an n by 3 matrix to polygonplot3d, were the rows of the matrix are the coordinates of the polygon's vertices.

In the third method we pass three vectors of length n each to polygonplot3d.  If these vectors are called U, V, W, then ( U[k], V[k], W[k] ) is taken to be the coordinates of the vertex number k.

plots:-polygonplot3d(U,V,W);

Your call to polygonplot3d invokes the thrid method. Thus, the first plotted point has coordinates [17/2,-17/2,0] as you have observed.  That's not what you want.

What you want is to invoke the first method which requires the vertices to be given as lists, not vectors.  There are two ways you may approach this: either provide vertices as lists, as in

A1 := [17/2, 0, 0];
A2 := [-17/2, 0, 0];
B1 := [0, 11/2, sqrt(166/4)];
plots:-polygonplot3d( [A1,A2,B1]);

or provide the vertices as vectors as you already have, and have maple convert them to lists before plotting, as in

A1 := Vector([17/2, 0, 0]);
A2 := Vector([-17/2, 0, 0]);
B1 := Vector([0, 11/2, sqrt(166/4)]);
plots:-polygonplot3d(convert~([A1,A2,B1], list));

For the weak formulation, you multiply the equation by a so-called "test function", let's say v, and then integrate the various terms by parts in order to reduce the differentiation order of the u function and increase the differentiation order of the v function.  For instance the fourth order derivative term, after integration by parts, changes to ∇² u . ∇²v .  I assume that you know how to do this, otherwise you will need to study some very basic ideas related to weak solutions of PDEs.  The formal computation above makes sense if the functions u and v are in the Sobolev space H². 

After you have done that, you introduce the Galerkin approximation of functions in H² whereby your PDE reduces to a system of linear ODEs in the time variable.

For the computational implementation you will need to pick basis functions for the Galerkin approximation.  That depends on the dimension of the space variable.  Is your space one or two or three dimensional? I am guessing that you are interested in the two-dimensional case.  In that case we will want triangular elements with Argyris shape functions.  (Do a search for "Argyris element" in google.)  Argyris shape functions are designed to connect across element boundaries as C¹ functions. That ensures that the Galerkin approximation forms a proper subspace of H².

• What is the unknown?  Is it u or w?
• What is f? Don't you mean f(x,y)?
• This being a second order equation in t, requires the specification of the initial velocity.
• The boundary conditions are insufficient.  For the biharmonic operator you need to specify not only the values of the unknown on the boundary, but also its derivatives.

Maple 2020 is unable to solve the following much simpler biharmonic equation in 1D:

restart;
pde := diff(w(x,t), t$2) + diff(w(x,t), x$4) = 0;
bc := w(0,t) = 0, w(1,t) = 0, D[1](w)(0,t)=0, D[1](w)(1,t)=0;
ic := w(x,0) = sin(Pi*x)^2, D[2](w)(x,0)=0;
pdsolve({pde,ic,bc});

I don't have Maple 2021. Give it a try with this 1D problem first.

PS: Maple 2020 can easily solve the 1D  problem above numerically. That doesn't help you with your 2D problem since Maple 2020's numeric PDE solver cannot handle PDEs with three independent variables.

The antiderivative of a function is determined up to an arbitrary additive constant.  So saying you want the antiderivative to be 0.739 is not meaningful.

Here is what you want:

restart;
f := sqrt(4+1/x)/2;
integ := int(f, x=0..a) assuming a > 0;
ans := fsolve(integ=0.739);
                      ans := 0.3726369981
# verify the answer
int(f, x=0..ans);
                          0.7390000000

As to your question about the shape of a chain draped over two pulleys, I suggest that you look first at the much simpler problem of determining the shape of a chain of length L suspended from two points within a vertical plane. Here is how it is done.

Let (x1,y1) and (x2,y2) be the Cartesian coordinates of the suspension points, and let y = f(x) be the equation of the chain.

The curve goes through the points (x1,y1) and (x2,y2).  That gives us two equations:
eq1 := y1 = f(x1),
eq2 := y2 = f(x2).

The curve's length is calculated by the usual curve length formula from calculus:
eq3 := int( sqrt(1 + f'(x)^2), x=x1..x2) = L;

Okay, so far we have specified the chain's geometry.  Now the physics.  It's possible to show, but I will not go through it here, that the chain's static equilibrium is expressed through the differential equation f''(x) = k*sqrt(1+f'(x)^2), where k is to be determined. In the worksheet below we see that this differential equation's general solution is a hyperbolic cosine.

Solving the differential equation brings along two integration constants c1 and c2.  So altogether, the differential equation's solution involves three unknown c1, c2, k.  We solve the system of equations {eq1, eq2, eq3} to determine those three unknowns.

The system is transcendental and has no useful symbolic solution, but we may solve the system numerically for any desired set of parameters (x1,y1), (x2,y2), and L.  See the worksheet below.

Going back to your question on the shape of a chain draped over two pulleys, that can be done in a similar way. The curve will be a catenary since that's the solution of the differential equation of the chain's equilibrium. I haven't done the calculations myself but the idea would be that instead of passing the curve from the points (x1,y1), (x2,y2), we need to find the curve which is tangent to the circles c1 and c2 that represent the pulleys. The solution should be within reach but not worth the trouble unless you really need it for some purpose.

Finally, regarding your question whether the points labeled A in your diagram in your worksheet are on horizontal diameters, the answer is no.  Just consider a chain which is tightly wrapped around the pulleys. The points A will be on vertical diameters in that case.

restart;
# The differential equation of a hanging chain
de := diff(y(x),x,x) = k*sqrt(1 + diff(y(x),x)^2);

# Solving the differential equation we obtain the equation of the hanging
# chain which involves three parameters, _C1, _C2, k:
dsolve(de):
f := unapply(rhs(%), x);

# We want the chain to go through the points (x1,y1) and (x2,y2):
eq1 := y1 = f(x1);
eq2 := y2 = f(x2);

# Additionally, we want the chain's length be a specified amount, L.
# The chain's length is calculated by the curve length formula from calculus:
int(sqrt(1+D(f)(x)^2), x):
simplify(%) assuming real:
eq3 := eval(%, x=x2) - eval(%, x=x1) = L;

# Now we need to solve the system of three equations {eq1, eq2, eq3} for the
# three unknowns _C1, _C2, k.  There is no symbolic solution, so we do  a
# numerical illustration.
params := { x1=0, y1=0, x2=10, y2=2, L=15 };
sys := eval({eq1,eq2,eq3}, params);
sol := fsolve(sys, {_C1,_C2,k}, {k=0..infinity});
#
# And now plot the solution
%plot(f(x), x=x1..x2, scaling=constrained):
eval(%, sol union params):
value(%);

Download worksheet: mw.mw

Here is the solution of your revised problem.  The solutions for f(x) and g(x) can be obtained entirely within Maple.

restart;

# Write the f derivative with the D operator
eq1 := x*D(f)(x^2) + diff(g(x),x) = cos(x) - 3*x^3;
eq2 := f(x^2) = sin(x) - x^4 - g(x);

# Calculate f'(x^2) from eq2 and substitute the result in eq1
# to obtain a differential equation in the unknown g(x):
diff(eq2, x);
isolate(%, D(f)(x^2));
de := subs(%, eq1);

# Solve the de for g(x), and then plug it in eq2 to calculate f(x);
dsolve(de);
eval(eq2, %);
algsubs(x^2=x, %);

The calculated f(x) and g(x) agree with what you have anticipated.

In your worksheet you refer to Problem #104 of 200 Puzzling Physics Problems. It is shown there that the shape of a chain in a steady motion between pulleys can be anything, so asking how to describe the chain's shape is pointless.

I must point out that the book's solution ignores gravity.  Adding gravity will drastically change the problem and make it quite nontrivial.  I don't expect a simple/elegant solution in that case.

The cantenary solution noted by Carl is good for a static chain.  As he has noted, a static uniform chain supported by two uprights of equal height assumes the shape of a simple hyperbolic cosine function. In fact, the "equal heights" assumption is unnecessary; the shape is a hyperbolic cosine even if the supports are at different levels.

Your candidates for f(x) and g(x) will work if you change x*f'(x) to x^2*f'(x) in the first equation:

restart;
eq1 := x^2*diff(f(x),x) + diff(g(x),x) = cos(x) - 3*x^3;
eq2 := f(x^2) + g(x) = sin(x) - x^4;
f := x -> -1/2*x^2 + C;
g := x -> sin(x) - 1/2*x^4 - C;

(lhs-rhs)(eq1);
                               0
(lhs-rhs)(eq2);
                               0

The above verifies that your candidates for f and g are correct.  I don't know how to solve the equations if f and g were not given. I don't expect that to be doable in Maple.

You have

pde33 := subs({D(f)(w)*(D^(3)(f)(w) = -c/3*(D^(4)(f)(w)}, pde3);

Those D^3 and D^4 should be D@@3 and D@@4.

After fixing a couple of misplaced parenthesis, the correct version of the above becomes

pde33 := subs({D(f)(w)*(D@@3)(f)(w) = -c/3*(D@@4)(f)(w)}, pde3);

The first version of the code that you posted is pretty much complete but it is carelessly written.  In Maple, there is a major difference between the "=" and ":=" operators.  The purpose of the first one is to state the equality of the left-hand and right-hand operands.  The second one is to define a symbol in terms of other symbols.   That's the first concept that you need to learn and master before you begin making use of Maple.

In the code that you have presented, those two operators are used randomly throughout, rendering the code meaningless.  For instance, where you have

x = a + i*h;

what you really mean is to define x to be a short-hand notation for a + i*h.  Therefore, that should be

x := a + i*h;

Your first task then is to go through that code very carefully and figure out where you need "=" and where ":=".

After you have fixed that, you need to change the few occurrences of rem(j,2).  What you need there is irem(j,2).  

If you do those two steps correctly, you will have a working proc.  Next, you need to learn how to use a proc.

Your proc is designed to integrate a function f(x,y) on a domain the x-y plane bounded by the lines x=a and x=b, and the curves y = c(x) and y=d(x).  To use the proc, you need to supply it with the three functions f, c, and d.  For instance, if f(x,y) = cos(y), c(x) = 0, and d(x) = x, you need

f := (x,y) -> cos(y);
c := x -> 0;
d := x -> x;
APV(f, 0.0, 1.0*Pi, c, d, 6, 6);

If you have been careful, this will produce the answer 2.000980454 which is pretty close to the exact solution which is 2.

I am not posting my corrected version of the code because I don't want to deprive you of the benefits of learning from this exercise.

You have α=0.1 and r goes from −1 to 1. What, in your opinion, is the value of r^(alpha−2) when r=−0.5?

I have fixed quite a few errors in your worksheet, most important of which are:

1. The ":=" and "=" mean quite different things to Maple.  The mathematical equal sign is "=".
2. The derivative at time t=0 should be written as D(xi)(0).

The numbers in your equations are quite odd.  Do you really want to deal with functions such as cos(10^12*t)?  That function oscillates billions of times in the time interval 0 < t < 1.

Anyway, here is the corrected worksheet, for whatever it's worth.
 

Download 66-corrected.mw