function x=gaussjordan(a,b)
  n=length(b);
  np1=n+1;
  carpan=0;
  C=a;
  for i=1:n
      C(i,n+1)=b(i);
  end
  for k=1:n
    %gauss jordan elimination
    for j=k+1:np1
      C(k,j)=C(k,j)/C(k,k);
    end
    for i=1:n
      if i~=k 
         carpan=C(i,k);
          C(i,k)=0;
          for j=k+1:np1
            C(i,j)=C(i,j)-carpan*C(k,j); end
       end    
    end
  end
  for i=1:n
   C(i,i)=1;   
   x(i)=C(i,np1);
  end
  end  

