function [B]=LU_band(A,m1,m2);
n=size(A,1);
B=zeros(n,m1+m2+1);
for i=m1:-1:1
    B(:,m1+1-i)=[zeros(i,1); diag(A,-i)];
end
B(:,m1+1)=diag(A);
for i=1:m2
    B(:,m1+1+i)=[diag(A,i);zeros(i,1)];
end

for i=1:n-1
    for j=(i+1):min(i+m1,n)
        B(j,i+m1+1-j)=B(j,i+m1+1-j)/B(i,m1+1);
        for k=(i+1):min(i+m2,n)
             B(j,k+m1+1-j)=B(j,k+m1+1-j)-B(j,i+m1+1-j)*B(i,k+m1+1-i);
        end
    end
end

%