function x = ThomasSolve(a, b, c, y)
% Funkcja rozwi¹zuje trójdiagonalny uk³ad równañ o M niewiadomych:
% [ b(1)  c(1)                                  ] [  x(1)  ]   [  y(1)  ]
% [ a(2)  b(2)  c(2)                            ] [  x(2)  ]   [  y(2)  ]
% [       a(3)  b(3)  c(3)                      ] [  x(3)  ]   [  y(3)  ]
% [            ...   ...   ...                  ] [  ...   ] = [  ...   ]
% [                    ...    ...    ...        ] [        ]   [        ]
% [                        a(M-1) b(M-1) c(M-1) ] [ x(M-1) ]   [ y(M-1) ]
% [                                 a(M)  b(M)  ] [  x(M)  ]   [  y(M)  ]
% gdzie y musi byæ wektorem (wierszowym) o d³ugoœci M, a, b oraz c musz¹
% byæ wektorami o d³ugoœci M (zauwa¿, ¿e a(1) i c(M) nie s¹ u¿ywane).

M = length(y);
x = zeros(M,1);
p = zeros(M,1);

w = b(1);
p(1) = -c(1)/w;
x(1) = y(1)/w;

for m = 2:M
    w = a(m)*p(m-1)+b(m);
    p(m) = -c(m)/w;
    x(m) = (y(m)-a(m)*x(m-1))/w;
end
for m = (M-1):-1:1
    x(m) = x(m) + p(m)*x(m+1);
end