function y=trid(a,b,c,d)
% solves tridiagonal system. a is the n-vector of diagonal elements.
% b subdiagonal, c superdiagonal are n-1 vectors, d right hand side.
n=length(a);
m=zeros(n,1);
y=zeros(n,1);
% factor matrix
f(1)=a(1);
for i=2:n
m(i)=b(i-1)/f(i-1);
f(i)=a(i)-m(i)*c(i-1);
end
% foward elimination
e(1)=d(1);
for i=2:n
e(i)=d(i)-m(i)*e(i-1);
end
% backsolve
y(n)=e(n)/f(n);
for k=n-1:-1:1
y(k)=(e(k)-c(k)*y(k+1))/f(k);
end
