function [u,w,l]=sls(Int,Bi,Bf,q,p,r)
%   
%   [u,w,l]=sls(Int,Bi,Bf,'q','p','r') uses finite differences to find
%   approximately the eigenvalues and eigen-functions of a Sturm-Liouville 
%   problem on the interval defined by the two-vector Int = [a,b] with 
%   coefficient functions defined in the m-files p.m, q.m, and r.m; i.e. 
%   we are solving the equation -(p*u')' + q*u = l*r*u for the function u 
%   and the eigenvalue l.  The initial and final homogeneous boundary 
%   conditions are specified in the-two vectors Bi and Bf.  For example if 
%   Bf=[c,d], then the boundary condition on the solution at the point b is 
%   c*u'(b) + d*u(b) = 0.  
%
%   The computation uses finite differences and sparse matrix routines.
%   The interval is divided into 1001 subintervals for the finite
%   difference approximation.  In the end the 20 smallest eigenvalues and
%   corresponiding eigenfunctions are determined.  The output consists 
%   of the matrix u with 1002 rows and 20 columns, each column consisting of 
%   the approximation to one of the eigenfunctions, the matrix w of the 
%   same size as u, the columns of which are the "dual" vectors to those of u, 
%   and the column vector l consisting of the approximations to the 
%   eigenvalues. 
%
%   The eigenfunctions are normalized by sum(u(:,j)^2) = 1 for each value of
%   j. 
%
%   The program allows the use of defaults as follows:
%   [u,w,l] = sls(Int,Bi,Bf,q,p) assumes r = 1.
%   [u,w,l] = sls(Int,Bi,Bf,q) assumes r = 1, and p = 1.
%   [u,w,l] = sls(Int,Bi,Bf) assumes r = 1, q = 0, and p = 1.
%
if (nargin < 3)
    error('There must be at least three input arguments.')
end
if (nargout < 3)
    error('There must be exactly three output arguments, e.g. [u,v,l] = sls(...).')
end
N = 1000;  % The number of interior division points.
K = 20;  % The number of eigenpairs computed.
a = Int(1);
b = Int(2);
h=(b-a)/(N+1);
x=a+[1:N]*h;
x2=a+h/2+[0:N]*h;
P = ones(1,N+1);
Q = zeros(1,N);
R = ones(1,N);
if (nargin > 5)
    P=feval(p,x2);
    Q=feval(q,x);
    R=feval(r,x);
    P(N+1)=feval(p,x2(N+1));
elseif (nargin == 5)
    P=feval(p,x2);
    Q=feval(q,x);
    P(N+1)=feval(p,x2(N+1));
    
elseif (nargin == 4)
    Q=feval(q,x);
    P(N+1)=1;
end


P1=P(1:N)';
P2=P(2:(N+1))';
P3=P(2:N)';
A = spdiags([[-P3;0],P1+P2+h^2*Q',[0;-P3]],-1:1,N,N);
% A = diag(P1 + P2) - diag(P3,1) - diag(P3,-1) + h*h*diag(Q);
if Bi(1)~=0
    if (nargin > 3) QQ=feval(q,a); else QQ = 0; end
    if (nargin > 4) PP=feval(p,a-h/2); else PP = 1; end
    if (nargin > 5) RR=feval(r,a); else RR = 1; end
    bb=(PP + P(1) + h*h*QQ - PP*2*h*Bi(2)/Bi(1))*P(1)/(PP+P(1));
    AA1=[bb,-P(1),zeros(1,N-1)];
    AA2=[-P(1),A(1,1:N)];
    AA3=[zeros(N-1,1),A(2:N,1:N)];
    A=[AA1;AA2;AA3];
    R=[RR*P(1)/(PP+P(1)),R];
end

if Bf(1)~=0
    if (nargin > 3) QQ=feval(q,a); else QQ = 0; end
    if (nargin > 4) PP=feval(p,a-h/2); else PP = 1; end
    if (nargin > 5) RR=feval(r,b); else RR = 1; end
    bb=(PP + P(N+1) + h*h*QQ + PP*2*h*Bf(2)/Bf(1))*P(N+1)/(PP+P(N+1));
    [c d]=size(A);
    AA1=[A(1:(c-1),1:c),zeros(c-1,1)];
    AA2=[A(c,:),-P(N+1)]; 
    AA3=[zeros(1,(c-1)),-P(N+1),bb];
    A=[AA1;AA2;AA3];
    R=[R,RR*P(N+1)/(PP+P(N+1))];
end

nr = length(R);
Rsp = spdiags(R',0,nr,nr);
B=h*h*Rsp;
A=B\A;
opts.disp = 0;
[V,D]=eigs(A,K,'sm',opts);
[l,k]=sort(spdiags(D));
u=V(:,k);
w=Rsp*u;
WW=diag(u'*w);
WW=sqrt(WW);
u=u/diag(WW);
w=w/diag(WW);

if Bi(1)==0
    [c d]=size(u);
    u=[zeros(1,d);u];
    w=[zeros(1,d);w];
end

if Bf(1)==0
    [c d] = size(u);
    u=[u;zeros(1,d)];
    w=[w;zeros(1,d)];
end