Legendre's polynomials and functions

I would like to improve this function to reduce the run time and optimize it.

function [Pl,Plm] = funlegendre(gamma)

    Plm = zeros(70,71);
    Pl = zeros(70,1);
    P0 = 1;
    Pl(1) = gamma;
    Plm(1,1) = sqrt(1-gamma^2);  

    for L=2:70
        for M=0:L

                if(M==0)
                    if (L==2)
                    Pl(L) = ((2*L-1)*gamma*Pl(L-1)-(L-1)*P0)/L;
                    else
                    Pl(l) = ((2*l-1)*gamma*Pl(l-1)-(l-1)*Pl(l-2))/l;
                    end
                elseif(M<L)
                    if(L==2)
                        if(M==1)
                        Plm(L,M) = (2*L-1)*Plm(1,1)*Pl(L-1);  
                        else
                        Plm(L,M) = (2*M-1)*Plm(1,1)*Plm(L-1,m-1); 
                        end
                    else
                        if(M==1)
                        Plm(L,M) = Plm(L-2,m) + (2*L-1)*Plm(1,1)*Pl(L-1);  
                        else
                        Plm(L,m) = Plm(L-2,m) + (2*L-1)*Plm(1,1)*Plm(L-1,M-1); 
                        end                 
                    end

                elseif(M==L)
                    Plm(L,M) = (2*L-1)*Plm(1,1)*Plm(L-1,L-1);
                else
                    Plm(L,M) = 0;
                end   
        end
    end
    Pl = sparse(Pl);
    Plm = sparse(Plm);
    end

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.

StackExchange Code Review Q#2743, answer score: 2