function [W1a,W1b,W2,W3]=eigbasU(m,N,x);

n=max(N,-N-1);
ind=1:length(x);
ind2=ind(abs(1-x.^2)>1e-6); %handle x=1 and x=-1 separately
indp=ind((abs(1-x.^2)<=1e-6)&(x>0));
indm=ind((abs(1-x.^2)<=1e-6)&(x<0));

    P=legendre(n,x);
    eta=1/sqrt(4*pi*gamma(n+m+1)/gamma(n-m+1)/(2*n+1));
    Ynm=(-1)^m*P(m+1,:);
    dYnm=zeros(size(Ynm));
    pYnm=zeros(size(Ynm));
    if(m==0)
        if(n==0) ; 
        else
            dYnm=P(2,:);
        end
    else
        dYnm(ind2)=-(-1)^m*(-m*x(ind2).*P(m+1,ind2)./sqrt(1-x(ind2).^2));
        if(m<n)
            dYnm(ind2)=dYnm(ind2)+(-1)^m*P(m+2,ind2);
        end
        pYnm(ind2)=Ynm(ind2)./sqrt(1-x(ind2).^2);
        if(m==1)
            dYnm(indp)=dYnm(indp)+n*(n+1)/2*ones(size(indp));
            dYnm(indm)=dYnm(indm)-(-1)^(n+1)*n*(n+1)/2*ones(size(indm));
            pYnm(indp)=dYnm(indp);
            pYnm(indm)=-dYnm(indm);
        end
        
    end
    
    % these are boundary values of the Lamb flows, normalized
    Unm1=eta*[Ynm; dYnm/(-N-1); sqrt(-1)*m*pYnm/(-N-1)]/sqrt((2*N+1)/(N+1));
    Unm3=eta*[Ynm; dYnm/N; sqrt(-1)*m*pYnm/N]/sqrt((2*N+1)/N);
    Unm2=eta*[zeros(size(Ynm)); sqrt(-1)*m*pYnm; -dYnm]/sqrt(N*(N+1));
    
   
    
    
% These are the corresponding Lamb flows
% W1a multiplies r^(N-1)
% W1b multiplies r^(N+1)
% W2 multiplies r^(N)
% W3 multiplies r^(N-1)
W1a=eta*[(2*N+3)/2*Ynm; (2*N+3)/(2*N)*dYnm; (2*N+3)/(2*N)*sqrt(-1)*m*pYnm]/sqrt((2*N+1)/(N+1));
W1b=eta*[-(2*N+1)/2*Ynm; -(2*N+1)*(N+3)/(2*N)/(N+1)*dYnm; 
    -(2*N+1)*(N+3)/(2*N)/(N+1)*sqrt(-1)*m*pYnm]/sqrt((2*N+1)/(N+1));
W2=Unm2;
W3=Unm3;

%check: W1a+W1b=Unm1