function [T,k,x,fa,d] = mysec(f,x0,x1,eps1,eps2)
%
% implementation of secant method
%

k=0
x = x0
fa = abs(f(x))
d=NaN
T(1,1:4)=[k,x,fa,d];
if fa < eps1
else
	k=1
    fa = abs(f(x1))
    d = abs(x1-x0)
    T(2,1:4)=[k,x1,fa,d];
	while fa >= eps1 & d >= eps2 & k <= 100
		k = k+1
		x= x1-f(x1)/(f(x1)-f(x0))*(x1-x0)
        x0 = x1
        x1 = x
	    fa=abs(f(x1))
		d=abs(x1-x0)
        T(k+1,1:4)=[k,x1,fa,d];
	end  
end

-----------------------------------------------------------------------------


function  [A,piv] = mylu(A)
% MYLU  LU factorization of an n-by-n matrix A.
%       mylu(A) returns M(1),...,M(n-1), P(1),...,P(n-1), and U such that
%       M(n-1)P(n-1)...M(1)P(1)A = U as described in Algorithm 3.4.1, 
p.111,
%       Golub and Van Loan(2nd Ed.). The output is stored as described 
there.
%       The upper triangle of A is overwritten with U while the spikes of
%       inverses of the M(i) are stored in the strict lower triangle of A.
%       The P(i) are represented by integers which are stored in
%       in a separate vector piv.

n = length(A);
piv = 1:n-1;
for k = 1:n-1
   [maxx,i] = max(abs(A(k:n,k)));
   if maxx ~= 0
      piv(k) = i + k - 1;
      A([k,piv(k)],1:n) = A([piv(k),k],1:n);
      A(k+1:n,k) = A(k+1:n,k)/A(k,k);
      A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
   end
end

------------------------------------------------------------------------------


function x = gaussjacobi(A,b,x0,eps)

n = length(A);
for i = 1:n
    C(i,:) = -A(i,:)/A(i,i);
    beta(i) = b(i)/A(i,i);
end
C = C + eye(n)
beta = beta'
x1 = C*x0+beta
k =1
while norm(x1-x0,inf)/norm(x1,inf) >= eps & k < 1000
    x0 = x1;
    x1 = C*x0 + beta
    k = k+1
end
x = x1


------------------------------------------------------------------------


function [x,C,beta] = gausseidel(A,b,x0,eps)

n = length(A);
for i = 1:n
    C(i,:) = -A(i,:)/A(i,i);
    beta(i) = b(i)/A(i,i);
end
C = C + eye(n)
beta = beta'
xh = x0;
for i = 1:n
    x1(i,1) = C(i,:)*xh+beta(i);
    xh(i,1) = x1(i,1);
    xh
end 
x1
k = 1
while norm(x1-x0,inf)/norm(x1,inf) >= eps && k < 1000
    x0 = x1;
    xh = x0;
    for i = 1:n
       x1(i,1) = C(i,:)*xh+beta(i);
       xh(i,1) = x1(i,1);
    end 
    x1
    k = k+1
end
x = x1