A = [ 1 1 1; 2 0 1; 0 -1 -1; -eye(3) ];
b = [ 15; 12; -5; 0; 0; 0];

% Add column to A with coefficients norm(A(i,:)) corresponding to variable
% r.
% A = [A, norms(A,2,2)];
Anorm = sqrt(sum(A.^2,2));
A = [A, Anorm];
[m n] = size(A);

% We want to maximize the radius, which is unknown x(4). Because linprog
% minimizes the objective, we use a negative cost coefficient.
c = [0;0;0;-1];
[x,fval,~,~,dual] = linprog(c,A,b);
y = dual.ineqlin;

% Display results.
fprintf(1,'The Chebyshev center coordinates are: \n');
disp(x(1:n-1));
fprintf(1,'The radius of the largest Euclidean ball is: \n');
disp(x(n));

% Which hyperplanes does the optimal ball touch?
ix = find(abs(A*x-b)<1e-6);
fprintf(1,'Hyperplanes touched by optimal ball\n')
disp(ix);

% The dual variables give us the change in optimal r for small changes
% in b.
fprintf(1,'Sensitivity of radius to small changes in b:\n')
disp(y);
