--- a/scripts/sparse/pcr.m
+++ b/scripts/sparse/pcr.m
@@ -110,7 +110,7 @@
 ## @sc{Example 1:} Simplest use of @code{pcr}
 ##
 ## @example
-##   x = pcr(A, b)
+##   x = pcr (A, b)
 ## @end example
 ##
 ## @sc{Example 2:} @code{pcr} with a function which computes
@@ -300,132 +300,143 @@
 
 endfunction
 
+
 %!demo
-%!
-%!      # Simplest usage of PCR (see also 'help pcr')
+%!  # Simplest usage of PCR (see also 'help pcr')
 %!
-%!      N = 20;
-%!      A = diag(linspace(-3.1,3,N)); b = rand(N,1); y = A\b; #y is the true solution
-%!      x = pcr(A,b);
-%!      printf('The solution relative error is %g\n', norm(x-y)/norm(y));
-%!
-%!      # You shouldn't be afraid if PCR issues some warning messages in this
-%!      # example: watch out in the second example, why it takes N iterations
-%!      # of PCR to converge to (a very accurate, by the way) solution
-%!demo
+%!  N = 20;
+%!  A = diag (linspace (-3.1,3,N)); b = rand (N,1);
+%!  y = A \ b;  # y is the true solution
+%!  x = pcr (A,b);
+%!  printf ("The solution relative error is %g\n", norm (x-y) / norm (y));
 %!
-%!      # Full output from PCR
-%!      # We use this output to plot the convergence history
+%!  # You shouldn't be afraid if PCR issues some warning messages in this
+%!  # example: watch out in the second example, why it takes N iterations
+%!  # of PCR to converge to (a very accurate, by the way) solution
+
+%!demo
+%!  # Full output from PCR
+%!  # We use this output to plot the convergence history
 %!
-%!      N = 20;
-%!      A = diag(linspace(-3.1,30,N)); b = rand(N,1); X = A\b; #X is the true solution
-%!      [x, flag, relres, iter, resvec] = pcr(A,b);
-%!      printf('The solution relative error is %g\n', norm(x-X)/norm(X));
-%!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)');
-%!      semilogy([0:iter],resvec/resvec(1),'o-g;relative residual;');
+%!  N = 20;
+%!  A = diag (linspace(-3.1,30,N)); b = rand (N,1);
+%!  X = A \ b;  # X is the true solution
+%!  [x, flag, relres, iter, resvec] = pcr (A,b);
+%!  printf ("The solution relative error is %g\n", norm (x-X) / norm (X));
+%!  clf;
+%!  title ("Convergence history");
+%!  xlabel ("Iteration"); ylabel ("log(||b-Ax||/||b||)");
+%!  semilogy ([0:iter], resvec/resvec(1), "o-g;relative residual;");
+
 %!demo
-%!
-%!      # Full output from PCR
-%!      # We use indefinite matrix based on the Hilbert matrix, with one
-%!      # strongly negative eigenvalue
-%!      # Hilbert matrix is extremely ill conditioned, so is ours,
-%!      # and that's why PCR WILL have problems
+%!  # Full output from PCR
+%!  # We use indefinite matrix based on the Hilbert matrix, with one
+%!  # strongly negative eigenvalue
+%!  # Hilbert matrix is extremely ill conditioned, so is ours,
+%!  # and that's why PCR WILL have problems
 %!
-%!      N = 10;
-%!      A = hilb(N); A(1,1)=-A(1,1); b = rand(N,1); X = A\b; #X is the true solution
-%!      printf('Condition number of A is   %g\n', cond(A));
-%!      [x, flag, relres, iter, resvec] = pcr(A,b,[],200);
-%!      if (flag == 3)
-%!        printf('PCR breakdown. System matrix is [close to] singular\n');
-%!      end
-%!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
-%!      semilogy([0:iter],resvec,'o-g;absolute residual;');
+%!  N = 10;
+%!  A = hilb (N); A(1,1) = -A(1,1); b = rand (N,1);
+%!  X = A \ b;  # X is the true solution
+%!  printf ("Condition number of A is   %g\n", cond (A));
+%!  [x, flag, relres, iter, resvec] = pcr (A,b,[],200);
+%!  if (flag == 3)
+%!    printf ("PCR breakdown.  System matrix is [close to] singular\n");
+%!  end
+%!  clf;
+%!  title ("Convergence history");
+%!  xlabel ("Iteration"); ylabel ("log(||b-Ax||)");
+%!  semilogy ([0:iter], resvec, "o-g;absolute residual;");
+
 %!demo
+%!  # Full output from PCR
+%!  # We use an indefinite matrix based on the 1-D Laplacian matrix for A,
+%!  # and here we have cond(A) = O(N^2)
+%!  # That's the reason we need some preconditioner; here we take
+%!  # a very simple and not powerful Jacobi preconditioner,
+%!  # which is the diagonal of A
 %!
-%!      # Full output from PCR
-%!      # We use an indefinite matrix based on the 1-D Laplacian matrix for A,
-%!      # and here we have cond(A) = O(N^2)
-%!      # That's the reason we need some preconditioner; here we take
-%!      # a very simple and not powerful Jacobi preconditioner,
-%!      # which is the diagonal of A
-%!
-%!      # Note that we use here indefinite preconditioners!
+%!  # Note that we use here indefinite preconditioners!
 %!
-%!      N = 100;
-%!      A = zeros(N,N);
-%!      for i=1:N-1 # form 1-D Laplacian matrix
-%!              A(i:i+1,i:i+1) = [2 -1; -1 2];
-%!      endfor
-%!      A = [A, zeros(size(A)); zeros(size(A)), -A];
-%!      b = rand(2*N,1); X = A\b; #X is the true solution
-%!      maxit = 80;
-%!      printf('System condition number is %g\n',cond(A));
-%!      # No preconditioner: the convergence is very slow!
+%!  N = 100;
+%!  A = zeros (N,N);
+%!  for i=1:N-1 # form 1-D Laplacian matrix
+%!    A(i:i+1,i:i+1) = [2 -1; -1 2];
+%!  endfor
+%!  A = [A, zeros(size(A)); zeros(size(A)), -A];
+%!  b = rand (2*N,1);
+%!  X = A \ b;  # X is the true solution
+%!  maxit = 80;
+%!  printf ("System condition number is %g\n", cond (A));
+%!  # No preconditioner: the convergence is very slow!
 %!
-%!      [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit);
-%!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
-%!      semilogy([0:iter],resvec,'o-g;NO preconditioning: absolute residual;');
-%!
-%!      pause(1);
-%!      # Test Jacobi preconditioner: it will not help much!!!
+%!  [x, flag, relres, iter, resvec] = pcr (A,b,[],maxit);
+%!  clf;
+%!  title ("Convergence history");
+%!  xlabel ("Iteration"); ylabel ("log(||b-Ax||)");
+%!  semilogy ([0:iter], resvec, "o-g;NO preconditioning: absolute residual;");
 %!
-%!      M = diag(diag(A)); # Jacobi preconditioner
-%!      [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
-%!      hold on;
-%!      semilogy([0:iter],resvec,'o-r;JACOBI preconditioner: absolute residual;');
+%!  pause (1);
+%!  # Test Jacobi preconditioner: it will not help much!!!
 %!
-%!      pause(1);
-%!      # Test nonoverlapping block Jacobi preconditioner: this one should give
-%!      # some convergence speedup!
+%!  M = diag (diag (A)); # Jacobi preconditioner
+%!  [x, flag, relres, iter, resvec] = pcr (A,b,[],maxit,M);
+%!  hold on;
+%!  semilogy ([0:iter],resvec,"o-r;JACOBI preconditioner: absolute residual;");
 %!
-%!      M = zeros(N,N);k=4;
-%!      for i=1:k:N # get k x k diagonal blocks of A
-%!              M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1);
-%!      endfor
-%!      M = [M, zeros(size(M)); zeros(size(M)), -M];
-%!      [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
-%!      semilogy([0:iter],resvec,'o-b;BLOCK JACOBI preconditioner: absolute residual;');
-%!      hold off;
+%!  pause (1);
+%!  # Test nonoverlapping block Jacobi preconditioner: this one should give
+%!  # some convergence speedup!
+%!
+%!  M = zeros (N,N); k = 4;
+%!  for i=1:k:N # get k x k diagonal blocks of A
+%!    M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1);
+%!  endfor
+%!  M = [M, zeros(size (M)); zeros(size(M)), -M];
+%!  [x, flag, relres, iter, resvec] = pcr (A,b,[],maxit,M);
+%!  semilogy ([0:iter], resvec, "o-b;BLOCK JACOBI preconditioner: absolute residual;");
+%!  hold off;
+
 %!test
-%!
-%!      #solve small indefinite diagonal system
-%!
-%!      N = 10;
-%!      A = diag(linspace(-10.1,10,N)); b = ones(N,1); X = A\b; #X is the true solution
-%!      [x, flag] = pcr(A,b,[],N+1);
-%!      assert(norm(x-X)/norm(X)<1e-10);
-%!      assert(flag,0);
+%! # solve small indefinite diagonal system
 %!
-%!test
-%!
-%!      #solve tridiagonal system, do not converge in default 20 iterations
-%!      #should perform max allowable default number of iterations
-%!
-%!      N = 100;
-%!      A = zeros(N,N);
-%!      for i=1:N-1 # form 1-D Laplacian matrix
-%!              A(i:i+1,i:i+1) = [2 -1; -1 2];
-%!      endfor
-%!      b = ones(N,1); X = A\b; #X is the true solution
-%!      [x, flag, relres, iter, resvec] = pcr(A,b,1e-12);
-%!      assert(flag,1);
-%!      assert(relres>0.6);
-%!      assert(iter,20);
-%!
+%! N = 10;
+%! A = diag (linspace (-10.1,10,N)); b = ones (N,1);
+%! X = A \ b;  # X is the true solution
+%! [x, flag] = pcr (A,b,[],N+1);
+%! assert (norm (x-X) / norm (X) < 1e-10);
+%! assert (flag, 0);
+
 %!test
-%!
-%!      #solve tridiagonal system with 'prefect' preconditioner
-%!      #converges in one iteration
+%! # solve tridiagonal system, do not converge in default 20 iterations
+%! # should perform max allowable default number of iterations
 %!
-%!      N = 100;
-%!      A = zeros(N,N);
-%!      for i=1:N-1 # form 1-D Laplacian matrix
-%!              A(i:i+1,i:i+1) = [2 -1; -1 2];
-%!      endfor
-%!      b = ones(N,1); X = A\b; #X is the true solution
-%!      [x, flag, relres, iter] = pcr(A,b,[],[],A,b);
-%!      assert(norm(x-X)/norm(X)<1e-6);
-%!      assert(relres<1e-6);
-%!      assert(flag,0);
-%!      assert(iter,1); #should converge in one iteration
+%! N = 100;
+%! A = zeros (N,N);
+%! for i=1:N-1 # form 1-D Laplacian matrix
+%!   A(i:i+1,i:i+1) = [2 -1; -1 2];
+%! endfor
+%! b = ones (N,1);
+%! X = A \ b;  # X is the true solution
+%! [x, flag, relres, iter, resvec] = pcr (A,b,1e-12);
+%! assert (flag, 1);
+%! assert (relres > 0.6);
+%! assert (iter, 20);
+
+%!test
+%! # solve tridiagonal system with "perfect" preconditioner
+%! # converges in one iteration
 %!
+%! N = 100;
+%! A = zeros (N,N);
+%! for i=1:N-1 # form 1-D Laplacian matrix
+%!   A(i:i+1,i:i+1) = [2 -1; -1 2];
+%! endfor
+%! b = ones (N,1);
+%! X = A \ b;  # X is the true solution
+%! [x, flag, relres, iter] = pcr (A,b,[],[],A,b);
+%! assert (norm (x-X) / norm(X) < 1e-6);
+%! assert (relres < 1e-6);
+%! assert (flag, 0);
+%! assert (iter, 1); # should converge in one iteration
+