Legendre's polynomials and functions
At first glance I can see a few points where you can start to optimize your polynom generation.
First thing I would refactor is the "heavy" use of conditional cases. There are only two well defined situation in which M==0 or M==L. Therefore you could extract both cases and implement them within the outer loop. If you don't care for code duplication copy/paste the functionality but if you care you have to write another function to call, which might be slow again. After you extracted the two special cases you cann ommit the remaining condition when you let L begin at 1 and run till L-1. Another point of optimization is to precompute values like 2*L-1 in the outer loop to have them "cached" for further access. As I'm not sure if sparse matrices are computationally faster than full matrices I'm not sure if I should recommend to use them from the beginning and not to convert the results into sparse matrices. Now I would try to vectorize as many operations as possible and let highly optimized code do the job for you.
And at last step @Elpezmuerto is fully right try to turn it into a mex function. This will precompile the function and speed up the execution time even further.