function x = tridiag(a,b,c,d)
% Solution of a system of linear equations with tridiagonal matrix
% We save just diagonals a, b, c:
% vector a = (0,   a_2,...,a_{n-1}, a_n) - lower diagonal
% vector b = (b_1, b_2,...,b_{n-1}, b_n) - main diagonal
% vector c = (c_1, c_2,...,c_{n-1}, 0)   - upper diagonal
% vector d = (d_1, d_2,...,d_{n-1}, d_n) - right hand side vector

% | b_1 c_1  0   
% | a_2 b_2 c_2  0   
% |  0  a_3 b_3 c_3  0
% |      0  a_4 b_4 c_4  0
% |          .   .   .   .

	n = length(d);
	x = zeros(n,1);

	c(1) = c(1)/b(1);
	d(1) = d(1)/b(1);
	b(1) = 1;

	for i=2:n
	    % elimination on  i-th row
	    b(i) = b(i)-c(i-1)*a(i);
	    d(i) = d(i)-d(i-1)*a(i);
	    a(i) = 0;
	    % division of the row by b(i)
	    c(i) = c(i)/b(i);
	    d(i) = d(i)/b(i);
	    b(i) = 1;
	end

	% back substitution:

	x(n) = d(n);
	for i=n-1:-1:1
	    x(i) = d(i)-c(i)*x(i+1);
	end 

end
