@dharr Thanks.
Sorry but since I'm very new with Maple and Groebner I did not succeed to do what you did on my 18 equations. Aren't my coefficients already rational? All my coefficients are with finate precision.
Anyway I tried multiplying all the diff_* equations by 10^8 in order to ensure integer coefficients and then called 'Groebner:-Solve(...)' for all my 18 equations. it was stuck and didn't return anything (had to stop the command).
Sorry for asking but can you to try to do what you did on my equations above?

Thanks

@vv 
Thanks for your answer.
I already solved a similar system of 10 polynomial equations (10 variables) and I used the 'solve' command that returned all the possible solutions. The 'solve' command returned at most 4 real solutions.
Also, when I solve this system (the 18 equations) with fsolve, I get many times the solution I expected to get after a few tries (and many times in the first try) so I believe that there are only a very few solutions to these equations (I would guess that it's never more than 10 solutions). so isn't there any trick that can be used to get all the solutions even though there are theoretically 2^18 solutions? maybe a smarter trick with the 'fsolve' command that will assure getting all the solutions?
BTW, the number 2^18 indeed sounds like a big number but how can I know if it's too big? (Maple was able to solve immediately a system of theoretically 2^10=1024 solutions). Do you know if there are other mathematical programs that maybe can solve it or you think it cannot be solved in any program?
Thanks

@acer 

I still don't understand why r22 can be any number between -1 to 1.

if I do all the substitutions that you said:

r13=0, r23=0, r31=0, r32=0, r33=1, t1=0, t2=0, t3=0, s1=0, s2=0, s3=0

then eq4 turns to be:

2*r22*sqrt(1-r22^2)=0

The solutions for this equation are r22=0, r22=1, r22=-1.

Meaning that the rotation matrices can be one of these 3 rotations:

R=Identity(3), R=[-1,0,0 ; 0,-1,0 ; 0,0,1], R=[0, 1, 0; 1, 0, 0 ; 0, 0, 1]

Where do you see that r22 can be any number between -1 to 1?

There is one additional constraint that I didn't include in the 15 equations.

That det(R) must be 1 (in the current equations det(R) can also be -1).

I meant to remove redundant solutions after I solve the equations. 

I tried to update my original question and to add the mathematical problem that I'm working on. I actually already published it here:

https://math.stackexchange.com/questions/2297050/finding-rotation-of-3-given-lines-in-3d-until-intersection-with-3-other-given-li

But here I tried to solve it in a different way using R*R^T=Identity (instead of dealing with trigonometric equations with 3 unknown angles).

But I'm not familiar with Mapleprimes and I didn't succeed to write the equations in latex form like I did in mathematica stackexchange.

@acer 

I'm not sure I understand why you gave values to some of the variables.

if you decided that r13=0, r23=0, r31=0, r32=0, r33=1, It means that the rotation matrix R=[r11 r12 r13; r21 r22 r23; r31 r32 r33] is actuall a rotation around the z axis (and also more possible matrices for the case that det(R)=-1). there are of course infinite such rotations. but I think that only the identity matrix will solve all the 15 equations in this case (meaning you should get that r11=1, r12=0, r21=0, r22=1).

Sorry I didn't post the mathematical question here (only the final equations that I got and I want Maple to solve). but the meaning of r_ij is actually elements in a rotation matrix.

Are you saying that there are infinite solutions? because I'm pretty sure there should be finite number of solutions to my original problem.

I didn't understand the last row of zeros (still new with Maple).

@one man 

as you said rotation matrix depends on 3 angles, therefore 3 degrees of freedom. therefore I have equations 1-6 (9 variables, 6 equations, meaning 3 degrees of freedom). The first 6 equations came from R*Transpose(R)=Identity_matrix (which is one of the properties of every rotaiton matrix). a rotation matrix must also satisfy that det(R)=1 but in my case I might get either det(R)=1 or det(R)=-1. I just meant to remove redundant solutions after I get all the solutions.

This is the definition of a polynomial equation of 3 variables with 2 degrees:

and I also wrote the definition in my original question:

"Meaning, every monomial has at most mulplications of 2 different variables or a single variable squared"

All my 15 equations follow that rull. So I still don't understand why you're saying that they are not polynomials.

About the mathematical problem that I'm trying to solve, I published it in Mathematica stack overflow. Here's the link if you're interested:

https://math.stackexchange.com/questions/2297050/finding-rotation-of-3-given-lines-in-3d-until-intersection-with-3-other-given-li

But here I tried to solve it in a different way (same problem exactly). I just came here to ask a question about Maple since I thought this is the goal of Mapleprimes and not to solve mathematical problems.

Anyway, I still will very glad to know if Maple is able to solve these kind of equations (and return ALL the solutions).

Thanks

@one man

The system of equations has exactly 15 variables. The 15 variables are:

r11,r12,r13,r21,r22,r23,r31,r32,r33,t1,t2,t3,s1,s2,s3.

You were looking at a specific example (of the solution) in which I gave values to the variables. in this example r22=r33, r23=-r32. but it doesn't have to be like that. it's not really the general values of r_ij. I just gave an example of a rotation matrix R=[r11 r12 r13;r21 r22 r23;r31 r32 r33] around the x axis. it can be also around y or z or any combination of them.

and why are you saying that they are not polynomials?

@tomleslie 

As I said in my original question, I know there's an algorithm for calculating ALL the solutions of any system of polynomial equations. But using 'fsolve' in Maple I always get only a single solution (a trivial solution. any solution in which all the values of t1,t2,t3,s1,s2,s3 are zeros is trivial and not good for me). and using 'solve' Maple is just not responding for hours (I'm not sure when to use 'solve' and when to use 'fsolve' since I always used 'solve' for simple numerical equations and it worked fine).

I want to know if I can some how get all the solutions to my 15 equations.

btw, I tried to thumb up your answer, and it didn't let me. It said that I need at least 10 reputations (even though I have exactly 10 reputations). I don't know why it didn't let me do that.

@tomleslie 

Thanks for your help. I didn't know that I should use fsolve rather than solve.

However, fsolve returns only a single solution (and not the solution I wrote in my original question that I expected to get).

Is there a way to get all the possible solutions to the equations?

@acer 

I wrote that excatly:

eq = maple('zip( `=`, exp1, exp2 )')

but in this command exp1, exp2 are considered new variables that I just defined (in maple) and not the variabales that I already defined in matlab.

I tried to do eq=eval(eq) but then matlab crushed.

The same with the 'children' command. I didn't succeed to make it work (in matlab) with maple engine even after you told me about the equivalent 'op' command in maple.

It's really frustrating that there's no link anywhere that explains how to work with this toolbox.