function [x_new,lambda_new,s_new]=quadip(n,eps)
% quadratic functional with inequality constraints solved with interior point method
% min 1/2 x^T A x -b^T x under the constraints x >= c
% the KKT conditions are (a slack variable s is added)
%  Ax-b-lambda=0
%  x-c-s=0
%  s_i lambda_i = 0 i=1,...,n
%  s>=0 and lambda>=0
% Newton method is used for the 3 first conditions
A = (n+1)*(n+1)*(sparse(2:n,1:n-1,-1,n,n) 
                + sparse(1:n,1:n,2,n,n) + sparse(1:n-1,2:n,-1,n,n));
b = -ones(n,1);
c = -0.05*ones(n,1);
x_old=zeros(n,1);
s_old=max(eps,ones(n,1));
lambda_old=max(eps,ones(n,1));
for iter=1:100
  mat1 = horzcat(???,-speye(n),sparse(n,n));
  mat2 = horzcat(speye(n),sparse(n,n),-speye(n));
  mat3 = horzcat(sparse(n,n),sparse(1:n,1:n,s_old,n,n),sparse(1:n,1:n,???,n,n));
  mat  = vertcat(mat1,mat2,mat3);
  rhs1 = A*x_old-b-???;
  rhs2 = x_old-c-???;
  rhs3 = sparse(n,1);
  for i=1:n
	  rhs3(i)=lambda_old(i)*???;
  end
  rhs = vertcat(rhs1,rhs2,rhs3);
  sol=mat\rhs;
  x_new=x_old-sol(1:n);
  plot(x_new);pause();
  lambda_new=max(eps,lambda_old-sol(n+1:???));
  s_new=max(eps,s_old-sol(???:3*n));
  discrep=norm(x_new-x_old)/norm(x_new);
  printf ("iter: %d Discrepancy: %f \n",iter,discrep);
  %fprintf("iter: %d Discrepancy: %f \n",iter,discrep); %matlab
  x_old=x_new;
  lambda_old=lambda_new;
  s_old=s_new;
  if (discrep<0.0001) 
    break
  end
end
end
