function [x] = solve3diag(d,l,u,b)
% function [x] = solve3diag(d,l,u,b)
%   Solve the system with 3-diagonal matrix A
%   d - main diagonal of A
%   l - subdiagonal of A
%   u - superdiagonal of A
%   b - the right side of the system
%   The variables d,l,u may be constant, it implies constant diagonals.
%   The variable b must be a vector because of the size of the system.
%

n=length(b);    % size of the system

% The test of sizes follows.
if length(d)==1, d=d*ones(1,n); end
if length(l)==1, l=l*ones(1,n-1); end
if length(u)==1, u=u*ones(1,n-1); end
if length(d)~=n || length(l)~=n-1 || length(u)~=n-1
    error('Incompatible size of input variables!');
end

x=zeros(n,1);

% Gaussian elimination follows
for k=1:n-1
    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 follows
x(n)=b(n)/d(n);
for k=n-1:-1:1
    x(k)=(b(k)-x(k+1)*u(k))/d(k);
end

end     % That's all, Folks!


