function table3_5

% tridaigonal
n=[4 8 12 16 160 1600 16000];
%n=[4 8 12 16 160];
N=length(n);
NN=length(n);
for in=1:N
    A=zeros(n(in),n(in));
    for i=1:n(in)
        A(i,i)=2;
        if i > 1
            A(i,i-1)=1;
        end
        if i < n(in)
            A(i,i+1)=1;
        end
    end
    xe=ones(n(in),1);
    a=2*ones(n(in),1); b=ones(n(in),1);  c=ones(n(in),1);
    f=4*ones(n(in),1); f(1)=3; f(n(in))=3;
    xc = tridiag( a, b, c, f );
    err=(xe-xc');
    eM(in)=norm(err,inf);
    cA(in)=cond(A,inf);
    ek(in)=eps*cond(A,inf);
end

% vandermonde
eMv=zeros(N); cAv=zeros(N);ekv=zeros(N);
%n=[4 8 12 16 160];
Nv=5;
for in=1:Nv
    A=zeros(n(in),n(in));
    for i=1:n(in)
        a=1/i;
        A(i,1)=1;
        for j=2:n(in)
            A(i,j)=a^(j-1);
        end
    end
    xe=ones(n(in),1);
    b=A*xe;
    xc=A\b;
    r=(b-A*xc);
    err=(xe-xc);
    eMv(in)=norm(err,inf);
    cAv(in)=cond(A,inf);
    ekv(in)=eps*cond(A,inf);
    %fprintf('\n %d & %17.2e & %17.1e & %17.1e',n(in),eM,cA,ek)
end

for in=1:NN
    fprintf('\n   %i  &  %5.1e  &  %5.1e  &  %5.1e      &  %5.1e  &  %5.1e  &  %5.1e  \\\\\\hline',n(in),eM(in),cA(in),ek(in),eMv(in),cAv(in),ekv(in))
end

fprintf('\n')


function y = tridiag( a, b, c, f )

%  Solve the  N x N  tridiagonal system for y:
%
%  [ a(1)  c(1)                                  ] [  y(1)  ]   [  f(1)  ]
%  [ b(2)  a(2)  c(2)                            ] [  y(2)  ]   [  f(2)  ]
%  [       b(3)  a(3)  c(3)                      ] [        ]   [        ]
%  [            ...   ...   ...                  ] [  ...   ] = [  ...   ]
%  [                    ...    ...    ...        ] [        ]   [        ]
%  [                        b(N-1) a(N-1) c(N-1) ] [ y(N-1) ]   [ f(N-1) ]
%  [                                 b(N)  a(N)  ] [  y(N)  ]   [  f(N)  ]
%
%  f must be a vector (row or column) of length N
%  a, b, c must be vectors of length N (note that b(1) and c(N) are not used)

N = length(f);
v = zeros(1,N);
y = v;

w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end

for j=N-1:-1:1
    y(j) = y(j) - v(j)*y(j+1);
end

