[nu1,nu2,mu1,mu2,theta,x,sx,id,pqa,ipqa,ierror]=dxpmu(nu1,nu2,mu1,mu2,theta,x,sx,id,pqa,ipqa,ierror);

function [nu1,nu2,mu1,mu2,theta,x,sx,id,pqa,ipqa,ierror]=dxpmu(nu1,nu2,mu1,mu2,theta,x,sx,id,pqa,ipqa,ierror);
persistent ip0 j mu n p0 x1 x2 ; 

if isempty(ip0), ip0=0; end;
if isempty(j), j=0; end;
if isempty(mu), mu=0; end;
if isempty(n), n=0; end;
%***BEGIN PROLOGUE  DXPMU
%***SUBSIDIARY
%***PURPOSE  To compute the values of Legendre functions for DXLEGF.
%            Method: backward mu-wise recurrence for P(-MU,NU,X) for
%            fixed nu to obtain P(-MU2,NU1,X), P(-(MU2-1),NU1,X), ...,
%            P(-MU1,NU1,X) and store in ascending mu order.
%***LIBRARY   SLATEC
%***CATEGORY  C3A2, C9
%***TYPE      doubleprecision (XPMU-S, DXPMU-D)
%***KEYWORDS  LEGENDRE FUNCTIONS
%***AUTHOR  Smith, John M., (NBS and George Mason University)
%***ROUTINES CALLED  DXADD, DXADJ, DXPQNU
%***REVISION HISTORY  (YYMMDD)
%   820728  DATE WRITTEN
%   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
%   901019  Revisions to prologue.  (DWL and WRB)
%   901106  Changed all specific intrinsics to generic.  (WRB)
%           Corrected order of sections in prologue and added TYPE
%           section.  (WRB)
%   920127  Revised PURPOSE section of prologue.  (DWL)
%***end PROLOGUE  DXPMU
if isempty(p0), p0=0; end;
if isempty(x1), x1=0; end;
if isempty(x2), x2=0; end;
pqa_shape=size(pqa);pqa=reshape(pqa,1,[]);
ipqa_shape=size(ipqa);ipqa=reshape(ipqa,1,[]);
%
%        CALL DXPQNU TO OBTAIN P(-MU2,NU,X)
%
%***FIRST EXECUTABLE STATEMENT  DXPMU
ierror = 0;
[nu1,nu2,mu2,theta,id,pqa,ipqa,ierror]=dxpqnu(nu1,nu2,mu2,theta,id,pqa,ipqa,ierror);
if( ierror==0 )
p0 = pqa(1);
ip0 = fix(ipqa(1));
mu = fix(mu2 - 1);
%
%        CALL DXPQNU TO OBTAIN P(-MU2-1,NU,X)
%
[nu1,nu2,mu,theta,id,pqa,ipqa,ierror]=dxpqnu(nu1,nu2,mu,theta,id,pqa,ipqa,ierror);
if( ierror==0 )
n = fix(mu2 - mu1 + 1);
pqa(n) = p0;
ipqa(n) = fix(ip0);
if( n~=1 )
pqa(n-1) = pqa(1);
ipqa(n-1) = fix(ipqa(1));
if( n~=2 )
j = fix(n - 2);
%
%        BACKWARD RECURRENCE IN MU TO OBTAIN
%              P(-MU2,NU1,X),P(-(MU2-1),NU1,X),....P(-MU1,NU1,X)
%              USING
%              (NU-MU)*(NU+MU+1.)*P(-(MU+1),NU,X)=
%                2.*MU*X*SQRT((1./(1.-X**2))*P(-MU,NU,X)-P(-(MU-1),NU,X)
%
while( true );
x1 = 2.0d0.*mu.*x.*sx.*pqa(j+1);
x2 = -(nu1-mu).*(nu1+mu+1.0d0).*pqa(j+2);
[x1,ipqa(j+1),x2,ipqa(j+2),pqa(j),ipqa(j),ierror]=dxadd(x1,ipqa(j+1),x2,ipqa(j+2),pqa(j),ipqa(j),ierror);
if( ierror~=0 )
break;
end;
[pqa(j),ipqa(j),ierror]=dxadj(pqa(j),ipqa(j),ierror);
if( ierror~=0 )
break;
end;
if( j==1 )
break;
end;
j = fix(j - 1);
mu = fix(mu - 1);
end;
end;
end;
end;
end;
pqa_shape=zeros(pqa_shape);pqa_shape(:)=pqa(1:numel(pqa_shape));pqa=pqa_shape;
ipqa_shape=zeros(ipqa_shape);ipqa_shape(:)=ipqa(1:numel(ipqa_shape));ipqa=ipqa_shape;
end %subroutine dxpmu
%DECK DXPMUP