startintlab; format long e; 

n = 20000; p_n = 224737; b = sparse(n,1); b(1) = 1;
A = spdiags(primes(p_n)',0,n,n); e = ones(n,2);
for k = 2.^(0:floor(log2(n))), A = A + spdiags(e,[-k k],n,n); end

% Session p. 161

tol = 1e-15; max_ite = 28;
D=spdiags(diag(A),0,n,n);
[x,fail,res,ite] = pcg(A,b,tol,max_ite,D);
if not(fail), x1 = midrad(x(1),norm(b-A*intval(x))), end

infsup(x1)