function [x error k] = BiCGSTAB(A, b, TOL, MAXITER)
	n=length(b); 
	TOL = TOL*norm(b); 
	error=[]; 
	x=zeros(n,1);
	rho=zeros(MAXITER+1,1);

	if norm(x)~=0
		r=b-A*x;
	else
		r=b;
	end
	
	hatr0=r;
	k=0; 
	rho(1)=1; 
	alpha=1; 
	omega=1;
	v=zeros(n,1);
	p=zeros(n,1);
	rho(2)=hatr0'*r;
	zeta=norm(r); 
	error=[error,zeta];

	while((zeta > TOL) & (k < MAXITER))
    		k=k+1;
    		if omega==0
       			error('Bi-CGSTAB breakdown, omega=0');
    		end
    		beta=(rho(k+1)/rho(k))*(alpha/omega);
    		p=r+beta*(p - omega*v);
    		v=A*p;
    		tau=hatr0'*v;
    		if tau==0
        			error('Bi-CGSTAB breakdown, tau=0');
    		end 
    		alpha=rho(k+1)/tau;
    		s=r-alpha*v; 
    		t=A*s;
    		tau=t'*t;
    		if tau==0
       			error('Bi-CGSTAB breakdown, t=0');
    		end
    		omega=t'*s/tau; 
    		rho(k+2) = -omega*(hatr0'*t);
%		rho(k+2) = hatr0'*r;
    		x=x+alpha*p+omega*s;
    		r=s-omega*t;
    		zeta=norm(r);
    		error=[error, zeta];
	end

end