function [x] = tridiag(a,b,c,f)
s = numel(a);
rho  = zeros(s+1,1);
mu  = zeros(s+1,1);
x  = zeros(s+1,1);
c(1) = 0;
b(s) = 0;
rho(1) = 0;
mu(1) = 0;
x(s+1) = 0;
for i=1:s
        mu(i+1) = -b(i)/(c(i)*mu(i)+a(i));
        rho(i+1)= (f(i)-c(i)*rho(i))/(c(i)*mu(i)+a(i));
end

for i = s:-1:1,
        x(i) = mu(i+1)*x(i+1) +rho(i+1);
end
x = x(1:s);
end