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

function [nu1,nu2,mu1,theta,x,sx,id,pqa,ipqa,ierror]=xqnu(nu1,nu2,mu1,theta,x,sx,id,pqa,ipqa,ierror);
persistent dmu ipq ipq1 ipq2 ipql1 ipql2 k mu nu pq pq1 pq2 pql1 pql2 x1 x2 ; 

if isempty(ipq), ipq=0; end;
if isempty(ipq1), ipq1=0; end;
if isempty(ipq2), ipq2=0; end;
if isempty(ipql1), ipql1=0; end;
if isempty(ipql2), ipql2=0; end;
if isempty(k), k=0; end;
if isempty(mu), mu=0; end;
%***BEGIN PROLOGUE  XQNU
%***SUBSIDIARY
%***PURPOSE  To compute the values of Legendre functions for XLEGF.
%            Method: backward nu-wise recurrence for Q(MU,NU,X) for
%            fixed mu to obtain Q(MU1,NU1,X), Q(MU1,NU1+1,X), ...,
%            Q(MU1,NU2,X).
%***LIBRARY   SLATEC
%***CATEGORY  C3A2, C9
%***TYPE      SINGLE PRECISION (XQNU-S, DXQNU-D)
%***KEYWORDS  LEGENDRE FUNCTIONS
%***AUTHOR  Smith, John M., (NBS and George Mason University)
%***ROUTINES CALLED  XADD, XADJ, XPQNU
%***REVISION HISTORY  (YYMMDD)
%   820728  DATE WRITTEN
%   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
%   901019  Revisions to prologue.  (DWL and WRB)
%   901106  Corrected order of sections in prologue and added TYPE
%           section.  (WRB)
%   920127  Revised PURPOSE section of prologue.  (DWL)
%***end PROLOGUE  XQNU
pqa_shape=size(pqa);pqa=reshape(pqa,1,[]);
ipqa_shape=size(ipqa);ipqa=reshape(ipqa,1,[]);
if isempty(dmu), dmu=0; end;
if isempty(nu), nu=0; end;
if isempty(pq), pq=0; end;
if isempty(pq1), pq1=0; end;
if isempty(pq2), pq2=0; end;
if isempty(x1), x1=0; end;
if isempty(x2), x2=0; end;
if isempty(pql1), pql1=0; end;
if isempty(pql2), pql2=0; end;
%***FIRST EXECUTABLE STATEMENT  XQNU
ierror = 0;
k = 0;
pq2 = 0.0;
ipq2 = 0;
pql2 = 0.0;
ipql2 = 0;
if( mu1~=1 )
mu = 0;
%
%        CALL XPQNU TO OBTAIN Q(0.,NU2,X) AND Q(0.,NU2-1,X)
%
[nu1,nu2,mu,theta,id,pqa,ipqa,ierror]=xpqnu(nu1,nu2,mu,theta,id,pqa,ipqa,ierror);
if( ierror~=0 )
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;
return;
end;
if( mu1==0 )
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;
return;
end;
k =fix((nu2-nu1+1.5));
pq2 = pqa(k);
ipq2 = fix(ipqa(k));
pql2 = pqa(k-1);
ipql2 = fix(ipqa(k-1));
end;
mu = 1;
%
%        CALL XPQNU TO OBTAIN Q(1.,NU2,X) AND Q(1.,NU2-1,X)
%
[nu1,nu2,mu,theta,id,pqa,ipqa,ierror]=xpqnu(nu1,nu2,mu,theta,id,pqa,ipqa,ierror);
if( ierror==0 )
if( mu1~=1 )
nu = nu2;
pq1 = pqa(k);
ipq1 = fix(ipqa(k));
pql1 = pqa(k-1);
ipql1 = fix(ipqa(k-1));
while( true );
mu = 1;
dmu = 1.;
%
%        FORWARD RECURRENCE IN MU TO OBTAIN Q(MU1,NU2,X) AND
%              Q(MU1,NU2-1,X) USING
%              Q(MU+1,NU,X)=-2.*MU*X*SQRT(1./(1.-X**2))*Q(MU,NU,X)
%                   -(NU+MU)*(NU-MU+1.)*Q(MU-1,NU,X)
%
%              FIRST FOR NU=NU2
%
while( true );
x1 = -2..*dmu.*x.*sx.*pq1;
x2 =(nu+dmu).*(nu-dmu+1.).*pq2;
[x1,ipq1,dumvar3,ipq2,pq,ipq,ierror]=xadd(x1,ipq1,-x2,ipq2,pq,ipq,ierror);
if( ierror~=0 )
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;
return;
end;
[pq,ipq,ierror]=xadj(pq,ipq,ierror);
if( ierror~=0 )
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;
return;
end;
pq2 = pq1;
ipq2 = fix(ipq1);
pq1 = pq;
ipq1 = fix(ipq);
mu = fix(mu + 1);
dmu = dmu + 1.;
if( mu>=mu1 )
break;
end;
end;
pqa(k) = pq;
ipqa(k) = fix(ipq);
if( k==1 )
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;
return;
end;
if( nu<nu2 )
break;
end;
%
%              THEN FOR NU=NU2-1
%
nu = nu - 1.;
pq2 = pql2;
ipq2 = fix(ipql2);
pq1 = pql1;
ipq1 = fix(ipql1);
k = fix(k - 1);
end;
%
%         BACKWARD RECURRENCE IN NU TO OBTAIN
%              Q(MU1,NU1,X),Q(MU1,NU1+1,X),....,Q(MU1,NU2,X)
%              USING
%              (NU-MU+1.)*Q(MU,NU+1,X)=
%                       (2.*NU+1.)*X*Q(MU,NU,X)-(NU+MU)*Q(MU,NU-1,X)
%
pq1 = pqa(k);
ipq1 = fix(ipqa(k));
pq2 = pqa(k+1);
ipq2 = fix(ipqa(k+1));
while( nu>nu1 );
k = fix(k - 1);
x1 =(2..*nu+1.).*x.*pq1./(nu+dmu);
x2 = -(nu-dmu+1.).*pq2./(nu+dmu);
[x1,ipq1,x2,ipq2,pq,ipq,ierror]=xadd(x1,ipq1,x2,ipq2,pq,ipq,ierror);
if( ierror~=0 )
break;
end;
[pq,ipq,ierror]=xadj(pq,ipq,ierror);
if( ierror~=0 )
break;
end;
pq2 = pq1;
ipq2 = fix(ipq1);
pq1 = pq;
ipq1 = fix(ipq);
pqa(k) = pq;
ipqa(k) = fix(ipq);
nu = nu - 1.;
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 xqnu
%DECK XRED