function [x,y] = solveBP_DIF(p,q,f,I,BV,n)
% function [x,y] = solveBP_DIF(p,q,f,I,BV,n)
%
% Solving the linear boundary problem by Finite Differences Method
% The equation: -(p(x)y')' + q(x)y = f(x)
% !!! Be careful of the sign '-' on the previous line. !!!
%
% p,q,f     - annonymous functions
% I=[a,b]   - interval
% BV        - boundary values, BV(1)=y(a), BV(2)=y(b)
% n         - number of subintervals, equally-spaced knots are generated
%
% x - knots
% y - approximate solution in the knots
%
% Finite Differences: 
% -(p(x)y'(x))' ~ -[p(x-h/2)y(x-h)-(p(x-h/2)+p(x+h/2))y(x)+p(x+h/2)y(x+h)]/h^2
%

ya=BV(1);yb=BV(2);  % boundary values
a=I(1);b=I(2);
h=(b-a)/n;
x=linspace(a,b,n+1);
x=x(:);                     % column
xx=x(2:n);                  % inner points
x1=(x(1:end-1)+x(2:end))/2; % midpoints of the subintervals
P=p(x1);
Q=q(xx);
F=f(xx);

d=(P(1:end-1)+P(2:end))/h^2+Q;  % main diagonal of the 3-diagonal system
l=-P(2:end-1)/h^2;              % subdiagonal
%u=-P(2:end-1)/h^2;             % superdiagonal
u=l;                            % superdiagonal is the same as l
b=F;                            % the right side of the system
b(1)=b(1)+P(1)*ya/h^2;          % y(a) is moved to the right side
b(end)=b(end)+P(end)*yb/h^2;    % y(b) is moved to the right side


% solving of the 3-diagonal system follows

yy=zeros(n-1,1);  % solution in inner points

% Gaussian elimination
for k=1:n-2
    q=l(k)/d(k);
    % l(k)=0;   % zeroing of the subdiagonal, not necessary
    d(k+1)=d(k+1)-q*u(k);
    b(k+1)=b(k+1)-q*b(k);
end

% back substitution
yy(n-1)=b(n-1)/d(n-1);
for k=n-2:-1:1
    yy(k)=(b(k)-yy(k+1)*u(k))/d(k);
end

y=[ya;yy;yb];

end     % That's all, Folks!

