It seems there is a negitive sign error here

SphericalY(3,-3,theta,phi)*conjugate(SphericalY(3,3,theta,phi)*sin(theta));
convert(%, LegendreP);
simplify(%) assuming real;

I think it has to do with any spherical harmonic of the form SphericalY(L,-L,theta,phi)

@Carl Love What does that mean? Is there something wrong with my code, or maple? 

Here's the plots

restart;
l:=1;
m:=-1;

plot3d(SphericalY(l,m,theta, phi)*conjugate(SphericalY(l,m,theta, phi)), phi=0..2*Pi,theta=0..Pi,coords=spherical);
#INCORRECT

restart;
Y:=(theta,phi)-> simplify(convert( SphericalY(l,m,theta, phi),elementary));
l:=1;
m:=-l;
plot3d(Y(theta,phi)*conjugate(Y(theta,phi)), phi=0..2*Pi,theta=0..Pi,coords=spherical);
#CORRECT

@rlopez But here I have lambda=+1, and mu=-1. So on the demoninator we have
(lambda-mu)! = (lambda - (-lambda))! = (2*lambda )! =2! and not a negitive integer. It seems the error always occurs when mu=-lambda (lambda being a positive integer)

Those are not the plots I expected. The spherical harmonics squared should display cylindrical symmetry. I did manage a solution though. 

restart;

Y:=(l,m,theta, phi)->simplify(convert( SphericalY(l,m,theta, phi),elementary));

for l from 1 to 2 do  
   m:=-l;  
   s:=(theta,phi) -> Y(l,m,theta, phi);
   plot3d(s(theta,phi)*conjugate(s(theta,phi));, phi=0..2*Pi,theta=0..Pi,coords=spherical);#sphereplot(f, phi = 0 .. 2*Pi, theta = 0 ..Pi,  style =     patchnogrid, scaling = constrained, axes = boxed);
end do;

But I do not understand why do I need to convert the SphericalY to an 'elementary' version of it? Yet this works, and the other did not...

thanks, that works perfectly!

thanks, that works perfectly!

@Carl Love If I do it symbolically yes, and I tried using 'map' and 'RootOf' but this didnt work, (it switch order of some values when plotting). 
 

also I will be making the matrix very large later, and its hard to solve the equation symbolically.

I have the same problem. Have you found a good solution yet?

Thanks!

I have the same problem. Have you found a good solution yet?

Thanks!

Thanks! thats a lot faster! 

Thanks! thats a lot faster! 

@acer Yes, in this case its 0. But in general it wont be. for example the value for SphericalY(2,0) and SphereicalY(0,0)... integrated w Vsi, is non-zero.

@acer Yes, in this case its 0. But in general it wont be. for example the value for SphericalY(2,0) and SphereicalY(0,0)... integrated w Vsi, is non-zero.

@Carl Love your right, as far as I can tell, the solutions may be so close together fsolve cant tell them apart. So
if I start at a solution of:
t=0.00005, theta=0.00005 and the REAL next solution is
t=0.0000500001, theta= 0.0000500001,

fsolve may skip it to find something like 
t=0.003, theta=0.003

(I call this teleporting) This only occurs when there exist more than one solution. In my case there are either 1 or 3 solutions which partially depend on the time (t). (Even if time allows for mulitple solutions, the trajectory of the ball may not)