function [x,y] = solveBP1(p,q,r,f,I,BV,n)
% function [x,y] = solveBP1(p,q,r,f,I,BV,n)
% Solving the linear boundary problem by Finite Differences Method
% The equation: p(x)y'' + q(x)y' + r(x)y=f(x)
% p,q,r,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:   y''(x) ~ [y(x-h)-2y(x)+y(x+h)]/h^2
%                        y'(x) ~ [y(x+h)-y(x-h)]/2h
%

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
P=p(xx);
Q=q(xx);
R=r(xx);
F=f(xx);

d=-2*P/h^2+R;                       % main diagonal of the 3-diagonal system
l=P(2:end)/h^2-Q(2:end)/(2*h);      % subdiagonal
u=P(1:end-1)/h^2+Q(1:end-1)/(2*h);  % superdiagonal
b=F;                                % the right side of the system
b(1)=b(1)-P(1)*ya/h^2+Q(1)*ya/(2*h);            % y(a) is moved to the right side
b(end)=b(end)-P(end)*yb/h^2-Q(end)*yb/(2*h);    % 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!

