[w,nroww,nrow,nbandl,nbandu,b]=bnslv(w,nroww,nrow,nbandl,nbandu,b);

function [w,nroww,nrow,nbandl,nbandu,b]=bnslv(w,nroww,nrow,nbandl,nbandu,b);
%***BEGIN PROLOGUE  BNSLV
%***SUBSIDIARY
%***PURPOSE  Subsidiary to BINT4 and BINTK
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (BNSLV-S, DBNSLV-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%  BNSLV is the BANSLV routine from
%        * A Practical Guide to Splines *  by C. de Boor
%
%  Companion routine to  BNFAC . It returns the solution  X  of the
%  linear system  A*X = B  in place of  B , given the LU-factorization
%  for  A  in the work array  W from BNFAC.
%
% *****  I N P U T  ******
%  W, NROWW,NROW,NBANDL,NBANDU.....Describe the LU-factorization of a
%        banded matrix  A  of order  NROW  as constructed in  BNFAC .
%        For details, see  BNFAC .
%  B.....Right side of the system to be solved .
%
% *****  O U T P U T  ******
%  B.....Contains the solution  X , of order  NROW .
%
% *****  M E T H O D  ******
%     (With  A = L*U, as stored in  W) the unit lower triangular system
%  L(U*X) = B  is solved for  Y = U*X, and  Y  stored in  B . Then the
%  upper triangular system  U*X = Y  is solved for  X  . The calcul-
%  ations are so arranged that the innermost loops stay within columns.
%
%***SEE ALSO  BINT4, BINTK
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   800901  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%***end PROLOGUE  BNSLV
%
persistent i j jmax middle nrowm1 ; 

if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(jmax), jmax=0; end;
if isempty(middle), middle=0; end;
if isempty(nrowm1), nrowm1=0; end;
w_shape=size(w);w=reshape([w(:).',zeros(1,ceil(numel(w)./prod([nroww])).*prod([nroww])-numel(w))],nroww,[]);
b_shape=size(b);b=reshape(b,1,[]);
%***FIRST EXECUTABLE STATEMENT  BNSLV
middle = fix(nbandu + 1);
if( nrow==1 )
b(1) = b(1)./w(middle,1);
else;
nrowm1 = fix(nrow - 1);
if( nbandl~=0 )
%                                 FORWARD PASS
%            FOR I=1,2,...,NROW-1, SUBTRACT  RIGHT SIDE(I)*(I-TH COLUMN
%            OF  L )  FROM RIGHT SIDE  (BELOW I-TH ROW) .
for i = 1 : nrowm1;
jmax = fix(min(nbandl,nrow-i));
for j = 1 : jmax;
b(i+j) = b(i+j) - b(i).*w(middle+j,i);
end; j = fix(jmax+1);
end; i = fix(nrowm1+1);
end;
%                                 BACKWARD PASS
%            FOR I=NROW,NROW-1,...,1, DIVIDE RIGHT SIDE(I) BY I-TH DIAG-
%            ONAL ENTRY OF  U, THEN SUBTRACT  RIGHT SIDE(I)*(I-TH COLUMN
%            OF  U)  FROM RIGHT SIDE  (ABOVE I-TH ROW).
if( nbandu>0 )
i = fix(nrow);
while (1);
b(i) = b(i)./w(middle,i);
jmax = fix(min(nbandu,i-1));
for j = 1 : jmax;
b(i-j) = b(i-j) - b(i).*w(middle-j,i);
end; j = fix(jmax+1);
i = fix(i - 1);
if( i<=1 )
break;
end;
end;
b(1) = b(1)./w(middle,1);
else;
%                                A  IS LOWER TRIANGULAR .
for i = 1 : nrow;
b(i) = b(i)./w(1,i);
end; i = fix(nrow+1);
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
return;
end;
end;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
end %subroutine bnslv
%DECK BQR