@jalal The command R[n,phi](u) rotates the vector u about the vector n through an angle phi.  It assumes that the axis of rotation goes through the origin.  In your worksheet your axis doesn't, and that's why you are getting an unexpected result.

Here I have translated the objects in your worksheet so that the rotation axis goes through the origin, and the result is a torus, as expected:

It is not too difficult to modify the proc R so that it performs a rotation about an arbitrarily positioned axis but perhaps the extra headache is not worth the trouble.  It is easier to set things up so that the rotation axis goes through the origin.

rotate-circle-animate_2.mw

@Kitonum There is no need for implicitplot3d for the second graph.  Just plot3d will do:

9*x^2*(3-z)^2+9*y*(3+z)^4 = 16*(-z^2+9)^2;
solve(%, y);
plot3d(%, x=-8..8, z=-3..3, view=-3..3);

or better yet:

9*x^2*(3-z)^2+9*y*(3+z)^4 = 16*(-z^2+9)^2;
solve(%, y);
plot3d([x,%,z], x=-8..8, z=-3..3, view=[-8..8, -3..3, -3..3]);

@vv Thank you very much for that construction.  It certainly produces the desired result.  I don't quite see how it works since I am not familiar with the Iterator package.  I will need to look into that.

@Kitonum Thanks.  That's very good.  I had looked at combinat:-permute, but didn't make the connection.

Hello @mmcdara I can't make sense of that worksheet.  Solving a nonlinear differential equation with finite elements is not exactly trivial.  I don't see anything that resembles a solution in there.

But as I wrote earlier, plugging the equation and boundary conditions into Maple produces a solution with no effort at all.

It looks to me that you have picked up a Maple code from somewhere and are attempting to apply it to a totally unrelated problem.  What are, for instance, CreateVector() and CreateMatrix()?  Furthermore, this code applies LinearSolve() to solve the resulting discretized system but your equations are nonlinear, so that won't work.

By the way, Maple's dsolve() solves your boundary value problem readily.  Why would you want to write a finite element code to solve it?

@nm Properly rounding 1.667 to two demimal places should produce 1.67 but yours produces 1.66:

fix_number(1.667,2);
        1.66

Perhaps that can be corrected with some small adjustment?

@ i see no worksheet.  Please try it again.

Upload your existing worksheet so that people can modify it as requested.

When you reply to this message, you will see a big fat green arrow in the toolbar.  Use that to upload the worksheet.

The message that you have posted, says:

If you experience any problems during installation, please visit the Maplesoft online technical support center at http://www.maplesoft.com/support.

Go to that web page and follow any of the given links there to contact technical support.

@Earl That's nice.  I converted your "Reply" to an "Answer" since it does offer an alternative approach.

@Carl Love Thanks for the compliment.  I had expected that a procedure for rotating a curve about an arbitrary axis would be included in some Maple library but as far as I can tell, it isn't.

@mmcdara Thanks for the kind words.

The calculation in my original worksheet extends almost without change to the general case of rotating a any curve (not necessarily a straight line) about a given line.  This is implemented as a procedure in the attached worksheet. 

Here is a sample obtained by rotating the helix <cos(t), sin(t), t/4> (shown in green) about the line <t, 0, t>. (shown in red). I have removed a part of the surface to make it possible to see the inside.

Download worksheet revolve-curve.mw

Although you have a system of three equations, you seem to be interested in plotting the system's vector field in the [r,phi] plane.  That is not possible in general. The vector field is 3-dimensional, so even when restricted to the the [r,phi] plane, the vectors will point outside that plane in general.  It's not a 2-dimensional picture.

@Vortex Unfortunately I know nothing about this area of physics.  I can only surmise from what you have said that f is a function of frequency, and the integration limits a and b are lower and upper bounds on the frequencies.  I have no idea what q represents.

It seems to me that the complete resolution of the question that you are asking depends on a full understanding of the physical phenomena that it models.  I am afraid I don't have that understanding.