%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Perform LU decomposition to tridiagonal matrices, and then solve linear     %
%   system by forward substitution and backward substitution                    %
%                                                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   The matrix and right-hand side
	f=exp(0:1:10); n=length(f);
    h=1;
    a=zeros(10,10);
    a(1,1)=2*2*h;
    for i=2:n-2
        a(i,i)=2*2*h;
        a(i,i-1)=h;
        a(i-1,i)=h;
    end
    a, pause
    for i=2:n-1
        b(i-1)=3*(f(i+1)-f(i))-3*(f(i)-f(i-1));
    end

%   Declare unknowns

    x=zeros(n-2,1); y=zeros(n-2,1);

%   LU factorization

    for k=1:(n-3)
        for i=(k+1):k+1
            a(i,k)=a(i,k)/a(k,k);
            for j=(k+1):k+1
                a(i,j)=a(i,j)-a(i,k)*a(k,j);
            end
        end
    end

%   Forward substitution

    y=b;
    for i=2:n-2
        for j=(i-1):(i-1)
            y(i)=y(i)-a(i,j)*y(j);
        end
    end

%   Backward substitution

    for i=(n-2):(-1):1
        for j=(i+1):min((i+1),n-2)
            y(i)=y(i)-a(i,j)*x(j);
        end
        x(i)=y(i)/a(i,i);
    end
    x