--- a/scripts/sparse/nonzeros.m
+++ b/scripts/sparse/nonzeros.m
@@ -7,11 +7,11 @@
 ## Returns a vector of the non-zero values of the sparse matrix @var{s}.
 ## @end deftypefn
 
-function t = nonzeros(s)
-  if issparse(s)
-    [i,j,t] = spfind(s);
+function t = nonzeros (s)
+  if (issparse (s))
+    [i, j, t] = spfind (s);
   else
-    [i,j,t] = find(s);
+    [i, j, t] = find (s);
   endif
 endfunction
 

--- a/scripts/sparse/normest.m
+++ b/scripts/sparse/normest.m
@@ -27,7 +27,7 @@
 ## @code{normest} to converge.
 ## @end deftypefn
 
-function [e1, c] = normest(A, tol)
+function [e1, c] = normest (A, tol)
   if (nargin < 2)
     tol = 1e-6;
   endif
@@ -35,14 +35,14 @@
     tol = eps
   endif
   if (ndims(A) != 2)
-    error("A must be a matrix");
+    error ("normest: A must be a matrix");
   endif 
-  maxA = max(max(abs(A)));
+  maxA = max (max (abs (A)));
   c = 0;
   if (maxA == 0)
     e1 = 0
   else
-    [m, n] = size(A);
+    [m, n] = size (A);
     B = A / maxA;
     Bt = B';
     if (m > n)
@@ -51,27 +51,27 @@
       Bt = tmp;
     endif
     e0 = 0;
-    x = randn(min(m,n),1);
-    e1 = norm(x);
+    x = randn (min (m, n), 1);
+    e1 = norm (x);
     x = x / e1;
-    e1 = sqrt(e1);
-    if (issparse(A))
-      while (abs(e1 - e0) > tol * e1)
+    e1 = sqrt (e1);
+    if (issparse (A))
+      while (abs (e1 - e0) > tol * e1)
 	e0 = e1;
 	x = B * (Bt * x);
-	e1 = norm(x);
+	e1 = norm (x);
 	x = x / e1;
-	e1 = sqrt(e1);
+	e1 = sqrt (e1);
 	c = c + 1;
       endwhile
     else
       B = B * Bt;
-      while (abs(e1 - e0) > tol * e1)
+      while (abs (e1 - e0) > tol * e1)
 	e0 = e1;
 	x = B * x;
-	e1 = norm(x);
+	e1 = norm (x);
 	x = x / e1;
-	e1 = sqrt(e1);
+	e1 = sqrt (e1);
 	c = c + 1;
       endwhile
     endif

--- a/scripts/sparse/pcg.m
+++ b/scripts/sparse/pcg.m
@@ -277,10 +277,10 @@
 	T = T(2:iter-2,2:iter-2);
 	l = eig(T);
 	eigest = [min(l), max(l)];
-	## fprintf (stderr, "PCG condest: %g\n", eigest(2)/eigest(1));
+	## fprintf (stderr, "pcg condest: %g\n", eigest(2)/eigest(1));
       else
 	eigest = [NaN, NaN];
-	warning ("PCG: eigenvalue estimate failed: iteration converged too fast.");
+	warning ("pcg: eigenvalue estimate failed: iteration converged too fast.");
       endif
     else
       eigest = [NaN, NaN];
@@ -308,19 +308,19 @@
   if (iter >= maxit-2)
     flag = 1;
     if (nargout < 2)
-      warning ("PCG: maximum number of iterations (%d) reached\n", iter);
-      warning ("The initial residual norm was reduced %g times.\n", 1.0/relres);
+      warning ("pcg: maximum number of iterations (%d) reached\n", iter);
+      warning ("the initial residual norm was reduced %g times.\n", 1.0/relres);
     endif
   elseif (nargout < 2)
-    fprintf (stderr, "PCG: converged in %d iterations. ", iter);
-    fprintf (stderr, "The initial residual norm was reduced %g times.\n",...
+    fprintf (stderr, "pcg: converged in %d iterations. ", iter);
+    fprintf (stderr, "the initial residual norm was reduced %g times.\n",...
 	     1.0/relres);
   endif
 
   if (! matrix_positive_definite)
     flag = 3;
     if (nargout < 2)
-      warning ("PCG: matrix not positive definite?\n");
+      warning ("pcg: matrix not positive definite?\n");
     endif
   endif
 endfunction

--- a/scripts/sparse/pcr.m
+++ b/scripts/sparse/pcr.m
@@ -279,20 +279,20 @@
   if (iter >= maxit-2)
     flag = 1;
     if (nargout < 2)
-      warning ("PCR: maximum number of iterations (%d) reached\n", iter);
-      warning ("The initial residual norm was reduced %g times.\n", 1.0/relres);
+      warning ("pcr: maximum number of iterations (%d) reached\n", iter);
+      warning ("the initial residual norm was reduced %g times.\n", 1.0/relres);
     endif
   elseif (nargout < 2 && ! breakdown)
-    fprintf (stderr, "PCR: converged in %d iterations. \n", iter);
-    fprintf (stderr, "The initial residual norm was reduced %g times.\n",
+    fprintf (stderr, "pcr: converged in %d iterations. \n", iter);
+    fprintf (stderr, "the initial residual norm was reduced %g times.\n",
 	     1.0/relres);
   endif
 
   if (breakdown)
     flag = 3;
     if (nargout < 2)
-      warning ("PCR: breakdown occured.\n");
-      warning ("System matrix singular or preconditioner indefinite?\n");
+      warning ("pcr: breakdown occured:\n");
+      warning ("system matrix singular or preconditioner indefinite?\n");
     endif
   endif
 

--- a/scripts/sparse/spconvert.m
+++ b/scripts/sparse/spconvert.m
@@ -27,18 +27,18 @@
 
 function s = spconvert (m)
 
-  if issparse(m)
+  if (issparse (m))
     s = m;
   else
-    sz = size(m);
-    if (nargin != 1 || !ismatrix(m) || !isreal(m) || length(sz) != 2 || 
-	(sz(2) != 3 && sz(2) != 4))
-      error ("spconvert: input matrix must be either sparse or a three or four column");
-      error ("           real matrix");
+    sz = size (m);
+    if (nargin != 1 || ! ismatrix (m) || ! isreal (m)
+	|| length (sz) != 2 || (sz(2) != 3 && sz(2) != 4))
+      error ("spconvert: argument must be sparse or real matrix with 3 or 4 columns");
     elseif (sz(2) == 3)
       s = sparse (m(:,1), m(:,2), m(:,3));
     else
       s = sparse (m(:,1), m(:,2), m(:,3) + 1i*m(:,4));
     endif
   endif
+
 endfunction

--- a/scripts/sparse/spdiags.m
+++ b/scripts/sparse/spdiags.m
@@ -53,33 +53,33 @@
 
 function [A, c] = spdiags (v, c, m, n)
 
-    if (nargin == 1 || nargin == 2)
-      ## extract nonzero diagonals of v into A,c
-      [i, j, v, nr, nc] = spfind (v);
-      if (nargin == 1)
-        c = unique(j-i);  # c contains the active diagonals
-      endif
-      ## FIXME: we can do this without a loop if we are clever
-      offset = max (min (c, nc-nr), 0);
-      A = zeros (min (nr, nc), length (c));
-      for k = 1:length(c)
-	idx = find (j-i == c(k));
-	A(j(idx)-offset(k),k) = v(idx);
-      end
-    elseif (nargin == 3)
-      ## Replace specific diagonals c of m with v,c
-      [nr,nc] = size (m);
-      B = spdiags (m, c);
-      A = m - spdiags (B, c, nr, nc) + spdiags (v, c, nr, nc);
-    else
-      ## Create new matrix of size mxn using v,c
-      [j, i, v] = find (v);
-      offset = max (min (c(:), n-m), 0);
-      j += offset(i);
-      i = j-c(:)(i);
-      idx = i>0 & i<=m & j>0 & j<=n;
-      A = sparse (i(idx), j(idx), v(idx), m, n);
-      
+  if (nargin == 1 || nargin == 2)
+    ## extract nonzero diagonals of v into A,c
+    [i, j, v, nr, nc] = spfind (v);
+    if (nargin == 1)
+      ## c contains the active diagonals
+      c = unique (j-i);
     endif
+    ## FIXME: we can do this without a loop if we are clever
+    offset = max (min (c, nc-nr), 0);
+    A = zeros (min (nr, nc), length (c));
+    for k = 1:length (c)
+      idx = find (j-i == c(k));
+      A(j(idx)-offset(k),k) = v(idx);
+    endfor
+  elseif (nargin == 3)
+    ## Replace specific diagonals c of m with v,c
+    [nr, nc] = size (m);
+    B = spdiags (m, c);
+    A = m - spdiags (B, c, nr, nc) + spdiags (v, c, nr, nc);
+  else
+    ## Create new matrix of size mxn using v,c
+    [j, i, v] = find (v);
+    offset = max (min (c(:), n-m), 0);
+    j += offset(i);
+    i = j-c(:)(i);
+    idx = i > 0 & i <= m & j > 0 & j <= n;
+    A = sparse (i(idx), j(idx), v(idx), m, n);
+  endif
 
 endfunction

--- a/scripts/sparse/speye.m
+++ b/scripts/sparse/speye.m
@@ -29,7 +29,7 @@
 ## is taken to be the size of the matrix to create.
 ## @end deftypefn
 
-function s = speye(m,n)
+function s = speye (m, n)
   if (nargin == 1)
     if (isvector (m) && length(m) == 2)
       n = m(2);
@@ -40,13 +40,13 @@
       error ("speye: invalid matrix dimension");
     endif
   else
-    if (!isscalar (m) || !isscalar (n))
+    if (! isscalar (m) || ! isscalar (n))
       error ("speye: invalid matrix dimension");
     endif
   endif
 
-  lo = min([m,n]);
-  s = sparse(1:lo,1:lo,1,m,n);
+  lo = min ([m, n]);
+  s = sparse (1:lo, 1:lo, 1, m, n);
 endfunction
 
 %!assert(issparse(speye(4)))

--- a/scripts/sparse/spfun.m
+++ b/scripts/sparse/spfun.m
@@ -23,22 +23,23 @@
 ## function handle or an inline function.
 ## @end deftypefn
 
-function t = spfun(f,s)
+function t = spfun (f, s)
+
   if (nargin != 2)
     print_usage ();
   endif
 
-  if issparse(s)
-    [i,j,v,m,n] = spfind(s);
+  if (issparse (s))
+    [i,j,v,m,n] = spfind (s);
   else
-    [i,j,v] = find(s);
-    [m,n] = size(s);
+    [i, j, v] = find (s);
+    [m, n] = size (s);
   end
 
   if (isa (f, "function_handle") || isa (f, "inline function"))
-    t = sparse(i,j,f(v),m,n);
+    t = sparse (i, j, f(v), m, n);
   else
-    t = sparse(i,j,feval(f,v),m,n);
+    t = sparse(i, j, feval (f, v), m, n);
   endif
 
 endfunction

--- a/scripts/sparse/sphcat.m
+++ b/scripts/sparse/sphcat.m
@@ -26,9 +26,9 @@
 
   persistent sphcat_warned = false;
 
-  if (!sphcat_warned)
+  if (! sphcat_warned)
     sphcat_warned = true;
-    warning ("sphcat: This function is depreciated. Use horzcat instead");
+    warning ("sphcat: This function is depreciated.  Use horzcat instead");
   endif
 
   y = horzcat (varargin{:});

--- a/scripts/sparse/spones.m
+++ b/scripts/sparse/spones.m
@@ -21,14 +21,14 @@
 ## sparse matrix with the same structure as @var{x}.
 ## @end deftypefn
 
-function s = spones(s)
-  if issparse(s)
-    [i,j,v,m,n] = spfind(s);
+function s = spones (s)
+  if (issparse (s))
+    [i, j, v, m, n] = spfind (s);
   else
-    [i,j,v] = find(s);
-    [m,n] = size(s);
+    [i, j, v] = find (s);
+    [m, n] = size (s);
   end
-  s = sparse(i,j,1,m,n);
+  s = sparse (i, j, 1, m, n);
 endfunction
 
 %!assert(issparse(spones([1,2;3,0])))

--- a/scripts/sparse/sprand.m
+++ b/scripts/sparse/sprand.m
@@ -30,26 +30,28 @@
 ##      2004-10-20      Texinfo help and copyright message 
 
 function S = sprand (m, n, d)
-  if nargin == 1
-    [i,j,v,nr,nc] = spfind(m);
-    S = sparse (i,j,rand(size(v)),nr,nc);
-  elseif nargin == 3
+  if (nargin == 1)
+    [i, j, v, nr, nc] = spfind (m);
+    S = sparse (i, j, rand (size (v)), nr, nc);
+  elseif (nargin == 3)
     mn = n*m;
-    k = round(d*mn); # how many entries in S would be satisfactory?
-    idx=unique(fix(rand(min(k*1.01,k+10),1)*mn))+1; 
-                # idx contains random numbers in [1,mn]
-  		# generate 1% or 10 more random values than necessary
-		# in order to reduce the probability that there are less than k
-		# distinct values;
-    		# maybe a better strategy could be used
-     		# but I don't think it's worth the price
-    k = min(length(idx),k);  # actual number of entries in S
-    j = floor((idx(1:k)-1)/m);
+    ## how many entries in S would be satisfactory?
+    k = round (d*mn);
+    idx = unique (fix (rand (min (k*1.01, k+10), 1) * mn)) + 1; 
+    ## idx contains random numbers in [1,mn]
+    ## generate 1% or 10 more random values than necessary in order to
+    ## reduce the probability that there are less than k distinct
+    ## values; maybe a better strategy could be used but I don't think
+    ## it's worth the price
+
+    ## actual number of entries in S
+    k = min (length (idx), k);
+    j = floor ((idx(1:k)-1)/m);
     i = idx(1:k) - j*m;
-    if isempty(i)
-      S = sparse(m,n);
+    if (isempty (i))
+      S = sparse (m, n);
     else
-      S = sparse(i,j+1,rand(k,1),m,n);
+      S = sparse (i, j+1, rand (k, 1), m, n);
     endif
   else
     print_usage ();

--- a/scripts/sparse/sprandn.m
+++ b/scripts/sparse/sprandn.m
@@ -21,27 +21,28 @@
 ## This program is public domain
 ## Author: Paul Kienzle <pkienzle@users.sf.net>
 
-function S = sprandn(m,n,d)
-  if nargin == 1
-    [i,j,v,nr,nc] = spfind(m);
-    S = sparse(i,j,randn(size(v)),nr,nc);
-  elseif nargin == 3
+function S = sprandn (m, n, d)
+  if (nargin == 1)
+    [i, j, v, nr, nc] = spfind (m);
+    S = sparse (i, j, randn (size (v)), nr, nc);
+  elseif (nargin == 3)
     mn = m*n;
-    k = round(d*mn);
-    idx=unique(fix(rand(min(k*1.01,k+10),1)*mn))+1; 
-                # idx contains random numbers in [1,mn]
-  		# generate 1% or 10 more random values than necessary
-		# in order to reduce the probability that there are less than k
-		# distinct values;
-    		# maybe a better strategy could be used
-     		# but I don't think it's worth the price
-    k = min(length(idx),k);  # actual number of entries in S
-    j = floor((idx(1:k)-1)/m);
+    k = round (d*mn);
+    idx = unique (fix (rand (min (k*1.01, k+10), 1) * mn)) + 1; 
+    ## idx contains random numbers in [1,mn]
+    ## generate 1% or 10 more random values than necessary in order to
+    ## reduce the probability that there are less than k distinct
+    ## values; maybe a better strategy could be used but I don't think
+    ## it's worth the price.
+
+    ## actual number of entries in S
+    k = min (length (idx), k);
+    j = floor ((idx(1:k)-1)/m);
     i = idx(1:k) - j*m;
-    if isempty(i)
-      S = sparse(m,n);
+    if (isempty (i))
+      S = sparse (m, n);
     else
-      S = sparse(i,j+1,randn(k,1),m,n);
+      S = sparse (i, j+1, randn (k, 1), m, n);
     endif
   else
     print_usage ();

--- a/scripts/sparse/sprandsym.m
+++ b/scripts/sparse/sprandsym.m
@@ -42,7 +42,7 @@
     n1 = m1 + rem (n, 2);
     mn1 = m1*n1;
     k1 = round (d*mn1);
-    idx1 = unique (fix(rand(min(k1*1.01,k1+10),1)*mn1))+1; 
+    idx1 = unique (fix (rand (min (k1*1.01, k1+10), 1) * mn1)) + 1; 
     ## idx contains random numbers in [1,mn] generate 1% or 10 more
     ## random values than necessary in order to reduce the probability
     ## that there are less than k distinct values; maybe a better
@@ -54,7 +54,7 @@
     n2 = ceil (n/2);
     nn2 = n2*n2;
     k2 = round (d*nn2);
-    idx2 = unique (fix (rand (min (k2*1.01, k1+10), 1)*nn2)) + 1; 
+    idx2 = unique (fix (rand (min (k2*1.01, k1+10), 1) * nn2)) + 1; 
     k2 = min (length (idx2), k2);
     j2 = floor ((idx2(1:k2)-1)/n2);
     i2 = idx2(1:k2) - j2*n2;
@@ -63,7 +63,7 @@
       S = sparse (n, n);
     else
       S1 = sparse (i1, j1+1, randn (k1, 1), m1, n1);
-      S = [tril(S1), sparse(m1,m1); ...
+      S = [tril(S1), sparse (m1, m1); ...
 	   sparse(i2,j2+1,randn(k2,1),n2,n2), triu(S1,1)'];
       S = S + tril (S, -1)';
     endif

--- a/scripts/sparse/spstats.m
+++ b/scripts/sparse/spstats.m
@@ -16,13 +16,14 @@
 ## @code{spstats (@var{s})} the zeros may disappear.
 ## @end deftypefn
 
-function [count,mean,var] = spstats(S,j)
-  if nargin < 1 || nargin > 2
+function [count, mean, var] = spstats (S, j)
+
+  if (nargin < 1 || nargin > 2)
     print_usage ();
   endif
 
-  if nargin == 1
-    [i,j,v] = find (S);
+  if (nargin == 1)
+    [i, j, v] = find (S);
   else
     v = S;    
     i = 1:length (v);
@@ -32,7 +33,7 @@
 
   count = spsum (sparse (i, j, 1, n, m));
   if (nargout > 1) 
-    mean = spsum(S) ./ count; 
+    mean = spsum (S) ./ count; 
   end
   if (nargout > 2) 
     ## FIXME Variance with count = 0 or 1?

--- a/scripts/sparse/spvcat.m
+++ b/scripts/sparse/spvcat.m
@@ -26,7 +26,7 @@
 
   persistent spvcat_warned = false;
 
-  if (!spvcat_warned)
+  if (! spvcat_warned)
     spvcat_warned = true;
     warning ("spvcat: This function is depreciated. Use vertcat instead");
   endif

--- a/scripts/sparse/treeplot.m
+++ b/scripts/sparse/treeplot.m
@@ -28,119 +28,149 @@
 function treeplot (Tree, NodeS, EdgeS)
 
   if (nargin < 1 || nargin > 3 || nargout > 0)
-    error ("treeplot: wrong number of input/output arguments");
+    print_usage ();
   else
-    if (!ismatrix(Tree) || size(Tree)(1) != 1 || !isnumeric(Tree) 
-        || !isvector(Tree) || any(Tree>length (Tree)))
+    if (! ismatrix (Tree) || rows (Tree) != 1 || ! isnumeric (Tree) 
+        || ! isvector (Tree) || any (Tree > length (Tree)))
       error ("treeplot: the first input argument must be a vector of predecessors");
-    else	
-    
-      NodeStyle = "0*;;";                  ## the inicialization of node end edge style
+    else
+      ## the inicialization of node end edge style
+      NodeStyle = "0*;;";
       EdgeStyle = "1;;";      
       if (nargin > 2)
         EdgeStyle = EdgeS;
         if (nargin > 1) 
-	  if ((length(findstr(NodeS,"*")) == 0) 
-	    && (length(findstr(NodeS,"+")) == 0) 
-	    && (length(findstr(NodeS,"x")) == 0) )
-	    NodeStyle = [NodeS,"o"];
+	  if (length (findstr (NodeS, "*")) == 0
+	      && length (findstr (NodeS, "+")) == 0
+	      && length (findstr (NodeS, "x")) == 0)
+	    NodeStyle = [NodeS, "o"];
 	  else
 	    NodeStyle = NodeS;
 	  endif
         endif
       endif
-      
+
       Tree = Tree(:)';		            ## make it a row vector
       NodNumber = length (Tree);            ## the count of nodes of the graph
-      ChildNumber = zeros (1,NodNumber+1);  ## the number of childrens
+      ChildNumber = zeros (1, NodNumber+1); ## the number of childrens
       
-      for i = 1:(NodNumber)       ## VecOfChild is helping vector which is used to speed up the
-                                  ##   choose of descendant nodes
+      for i = 1:NodNumber
+        ## VecOfChild is helping vector which is used to speed up the
+        ## choose of descendant nodes
 
         ChildNumber(Tree(i)+1) = ChildNumber(Tree(i)+1) + 1;
       endfor
       Pos = 1;
-      for i = 1:(NodNumber+1)
+      for i = 1:NodNumber+1
         Start(i) = Pos;
 	Help(i) = Pos;
-	Pos = Pos + ChildNumber(i);
+	Pos += ChildNumber(i);
 	Stop(i) = Pos;
       endfor
       for i = 1:NodNumber        
         VecOfChild(Help(Tree(i)+1)) = i;  
-	Help(Tree(i)+1)=Help(Tree(i)+1)+1;
-      endfor                    ## VecOfChild is helping vector which is used to speed up the
-                                ##   choose of descendant nodes
-				
-      ParNumber = 0;            ## the number of "parent" (actual) node (it's descendants will 
-                                ##   be browse in the next iteration)
-      LeftMost = 0;		## the x-coordinate of the left most descendant of "parent node"
-                                ##   this value is increased in each leaf		
-      Level = NodNumber;        ## the level of "parent" node (root level is NodNumber)
-      Max = NodNumber;          ## NodNumber - Max is the height of this graph
-      St = [-1,0];		## main stack - each item consists of two numbers - the number of node and
-                                ##   the number it's of parent node
-				##   on the top of stack there is "parent node"
-      Skelet = 0;		## stack which is use to draw the graph edge (it 
-                                ##   have to be uninterupted line)
-      while (ParNumber!=-1)     ## the top of the stack
+	Help(Tree(i)+1) = Help(Tree(i)+1)+1;
+      endfor
+
+      ## the number of "parent" (actual) node (it's descendants will be
+      ## browse in the next iteration)
+      ParNumber = 0;
+
+      ## the x-coordinate of the left most descendant of "parent node"
+      ## this value is increased in each leaf		
+      LeftMost = 0;
+
+      ## the level of "parent" node (root level is NodNumber)
+      Level = NodNumber;
+
+      ## NodNumber - Max is the height of this graph
+      Max = NodNumber;
+
+      ## main stack - each item consists of two numbers - the number of
+      ## node and the number it's of parent node on the top of stack
+      ## there is "parent node"
+      St = [-1,0];
+
+      ## stack which is use to draw the graph edge (it have to be
+      ## uninterupted line)
+      Skelet = 0;
+
+      ## the top of the stack
+      while (ParNumber != -1)
 	if (Start(ParNumber+1) < Stop(ParNumber+1))
 	  idx = VecOfChild(Start(ParNumber+1):Stop(ParNumber+1)-1);
 	else
-	  idx = zeros(1,0);
-	endif                   ## add to idx the vector of parent descendants
+	  idx = zeros (1, 0);
+	endif
+        ## add to idx the vector of parent descendants
 	St = [St ; [idx', ones(fliplr(size(idx)))*ParNumber]];
-	                        ## add to stack the records relevant to parent descandant s
+	## add to stack the records relevant to parent descandant s
 	if (ParNumber != 0)
-	  Skelet = [Skelet ; ([ones(size(idx))*ParNumber; idx])(:)];
+	  Skelet = [Skelet; ([ones(size(idx))*ParNumber; idx])(:)];
 	endif
-	if (St(end,2) != ParNumber)  ## if there is not any descendant of "parent node":
+
+	## if there is not any descendant of "parent node":
+	if (St(end,2) != ParNumber)
 	  LeftMost = LeftMost + 1;
           XCoordinateR(ParNumber) = LeftMost;           
 	  Max = min (Max, Level);
-          if ((length(St)>1) && (find((shift(St,1)-St) == 0) >1) && St(end,2) != St(end-1,2))
-	                             ## return to the nearest branching
-            Position = (find((shift(St(:,2),1)-St(:,2)) == 0))(end)+1;      ## the position to return
-                          ## position is the position on the stack, where should be started 
-			  ## further search (there are two nodes which has the same parent node)
+          if ((length(St)>1) && (find((shift(St,1)-St) == 0) >1)
+	      && St(end,2) != St(end-1,2))
+	    ## return to the nearest branching the position to return
+	    ## position is the position on the stack, where should be
+	    ## started further search (there are two nodes which has the
+	    ## same parent node)
+            Position = (find((shift(St(:,2),1)-St(:,2)) == 0))(end)+1;
             ParNumberVec = St(Position:end,2);
-                          ## the vector of removed nodes (the content of stack form position to end)
+            ## the vector of removed nodes (the content of stack form
+	    ## position to end)
             Skelet = [Skelet; flipud(ParNumberVec)];
             Level = Level + length(ParNumberVec);
-	                  ## the level have to be decreased
+	    ## the level have to be decreased
             XCoordinateR(ParNumberVec) = LeftMost;
             St(Position:end,:) = [];
           endif	
-	  St(end,:) = [];       	## remove the next node from "searched branch"
-          ParNumber = St(end,1);	## choose new "parent node"
-	  if (ParNumber != -1)		## if there is another branch start to search it
+       	  ## remove the next node from "searched branch"
+	  St(end,:) = [];
+	  ## choose new "parent node"
+          ParNumber = St(end,1);
+	  ## if there is another branch start to search it
+	  if (ParNumber != -1)
 	    Skelet = [Skelet ; St(end,2); ParNumber];
             YCoordinate(ParNumber) = Level;	
 	    XCoordinateL(ParNumber) = LeftMost + 1;
 	  endif
-	else                      ## there were descendants of "parent nod"
-          Level = Level - 1;      ## choose the last of them and go on through it
+	else
+          ## there were descendants of "parent nod" choose the last of
+	  ## them and go on through it
+          Level--;
 	  ParNumber = St(end,1);
 	  YCoordinate(ParNumber) = Level;     
 	  XCoordinateL(ParNumber) = LeftMost+1;
 	endif
       endwhile
 
-      XCoordinate = (XCoordinateL + XCoordinateR) / 2;     ## calculate the x coordinates 
-               ## (the known values are the position of most left and most right descendants)
+      ## calculate the x coordinates (the known values are the position
+      ## of most left and most right descendants)
+      XCoordinate = (XCoordinateL + XCoordinateR) / 2;
 
-      axis ([0.5 LeftMost+0.5 Max-0.5 NodNumber-0.5], "nolabel");  ## set axis and graph size 
-      
-      plot (XCoordinate,YCoordinate,NodeStyle);           ## plot grah nodes 
       hold ("on");
+
+      ## plot grah nodes 
+      plot (XCoordinate,YCoordinate,NodeStyle);
       
-      Skelet = [Skelet; 0];     ## helping command - usable for plotting edges
+      ## helping command - usable for plotting edges
+      Skelet = [Skelet; 0];
       
-      idx = find (Skelet == 0);                           ## draw graph edges 
+      ## draw graph edges 
+      idx = find (Skelet == 0);
        
-      for i = 2:length(idx)                               ## plot each tree component in one loop
-	istart = idx(i-1) + 1;                            ## tree component start
-	istop = idx(i) - 1;                               ## tree component end
+      ## plot each tree component in one loop
+      for i = 2:length(idx)
+        ## tree component start
+	istart = idx(i-1) + 1;
+        ## tree component end
+	istop = idx(i) - 1;
 	if (istop - istart < 1)                          
 	  continue;
 	endif
@@ -148,6 +178,9 @@
 	      YCoordinate(Skelet(istart:istop)), EdgeStyle)
       endfor
       
+      ## set axis and graph size 
+      axis ([0.5, LeftMost+0.5, Max-0.5, NodNumber-0.5], "nolabel");
+      
       hold ("off");
       
     endif