Is there a reason for looking at xhat as a function of all those variables?  Why not take xhat to be a function of a_1 only.  The solution will involve the remaining variables as parameters.

Download mw.mw

There is no meaningful way to do what you are asking, and that's not because of any shortcoming on the part of Maple.

Your equations are defined in a three-dimensional phase space, that is, the associated vector field is 3D and the solution curve is a curve in 3D.  You plot the projection of that 3D curve onto the (V,In) plane.  That works fine since the projection of a 3D curve onto a 2D plane is still a curve.

In contrast, the projection of a 3D vector field onto a 2D plane is not a vector field because each point of your 2D plotting domain is the foot of the projection of all the vectors in 3D that lie "above" that point.

Change
if not is({args}, set(numeric)) then return ('procname')('args') end if;
to
if not type({args}, set(numeric)) then return ('procname')('args') end if;
and it will work.

Comment added later:

It is likely that numeric is not what you want.   For instance, Pi is not numeric.  To allow for symbolic constants such as Pi as arguments, change numeric to realcons.  Here is the way I would do it:

if not type([args], list(realcons)) then return 'procname'(args) end if;

I have changed your notation so that the variables are now named x[1], x[2], etc., and the equations are dxdt[1], dxdt[2], etc.  This makes it easy to address them. Working with 2D input makes me mad, so I have converted your document to 1D math.

Download mw.mw

Download mw.mw

This website does not allow me to post the contents of a Maple worksheet (something's broken?) so here is the text version of what I wanted to post.

We look at a nonlinear boundary value problem which has infinitely many solutions.  We specify three initial guesses which pick three different solutions.

restart;
de := diff(u(x),x,x) + surd(u(x),3)= 0;
bc := u(0) = 0, u(1) = 0;

# Solution #1
dsol := dsolve({de,bc}, numeric, method=bvp[midrich], approxsoln=[u(x)=x*(1-x)], maxmesh=200);
plots:-odeplot(dsol);

# Solution #2
dsol := dsolve({de,bc}, numeric, method=bvp[midrich], approxsoln=[u(x)=-x*(1-x)], maxmesh=200);
plots:-odeplot(dsol);

# Solution #3
dsol := dsolve({de,bc}, numeric, method=bvp[midrich], approxsoln=[u(x)=x*(1-x)*(x-1/2)], maxmesh=500);
plots:-odeplot(dsol);

[inserted by moderator]

Download _BVP_with_infnitely_many_solutions.mw

restart;
with(plots):
f := (x,y) -> x^4 - 3*x^2 - 2*y^3 + 3*y + 0.5*x*y; 
solve( {diff(f(x,y),x)=0, diff(f(x,y),y)=0}, {x,y}):
(p->eval([x,y,f(x,y)],p))~([%]);
pointplot3d(%, symbol=solidsphere, symbolsize=20, color=red):
plot3d(f(x,y), x=-2..2, y=-2..2, style=surfacecontour, contours=20):
display(%%,%, view=-5..2, lightmodel=light1, orientation=[-90,0,0]);

Your worksheet is suffering from the 2D-math-input-syndrome.  I waded through it (painfully) and corrected a few errors here and there.  Here is the result.  Normally I wouldn't touch a 2D document with a 10-foot pole.

Download dsol.mw

You must realize that what you have posted is not a good way of sharing your work with others.  Why not post a worksheet file (*.mw) like what most other people do?

And if you really want helpful input, you should document your work.  Explain what every step of the calculation is meant to do. Otherwise saying "it doesn't work" makes no sense since we don't know what you have in your head.

That said, I attempted to make sense of the mess that you have posted. It appears to me that the reason for "it doesn't work" is that at the very end you have

display([Head(0.3), Rod(), DikA, DikB, Crankpin(alpha)], ...

I think that the 0.3 in Head(0.3) should be a variable, but since you have not explained your notation, I have no idea what it is supposed to be.

I gather from what I could understand, that you want to do something like this:

See how that's done in this worksheet:  mw.mw

This does exactly what you want.

restart;
with(Physics[Vectors]):
Typesetting:-Settings(typesetdot=true):
r_(t) = ChangeBasis(rhs(r_(t) = x(t)*_i + y(t)*_j), cylindrical, alsocomponents);
OM_ := r(t)*_rho(t);
v_ := diff(OM_, t);
a_ := diff(v_, t);

Added later: Actually the line r_(t) = Change... is not needed.  Delete.

The number that you have calculated for the answer to part (a) is correct but I cannot see the logic behind that calculation.  You may want to explain your reasoning.  Here is my calculation:

params := rho=1000, g=9.81, R=1/2, H=2;   # H=height, R=radius
W := int((H-y)*rho*g*Pi*(R/H*y)^2, y=0..H);
eval(W, {params});

which yields 2568.251995, the same as yours.  For parts (b) and (c) I get 13080π and 39240π.

Printing that is the easy part:

What you want to do with the result is something else.

@erik10 See if this is what you want.

Note: The message on the output says "number of chars = 22".  The "22" is shown properly in the worksheet but comes out broken on this website. Not my fault.

Download mw.mw

Download mw.mw

If you solve the equation x^2 + y^2 = 4 for y, you will get two solutions, y=±sqrt(4 - x^2). The plus sign gives the equation of the upper half of the circle, and the minus sign gives the equation of the lower half. So the answer to Part (a) is y = sqrt(4 - x^2), and the plot is produced by

plot(sqrt(4 - x^2), x=-2..2, scaling=constrained);

As to part (b): Rotating that semicircle about the x axis produces a sphere of radius 2.

As to part (c): The volume enclosed by the sphere is given by the formula V_x that you have shown in your post. We have f(x) = sqrt(4 - x^2), and therefore f(x)^2 = 4 - x^2.  So the volume is:

V__x := Pi*int( 4 - x^2, x = -2 .. 2);