function x=thomas(A,b)
%risolve il sistema Ax=b con l'alg. di Thomas
%A dovrebbe essere tridiagonale

n=length(b);
x=zeros(n,1);
for i=2:n
    m=-A(i,i-1)/A(i-1,i-1);
    A(i,i-1)=0;
    A(i,i)=A(i,i)+m*A(i-1,i);
    b(i)=b(i)+m*b(i-1);
end
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
    x(i)=(b(i)-A(i,i+1)*x(i+1))/A(i,i);
end
