new file mode 100644
--- /dev/null
+++ b/doc/interpreter/sparseimages.m
@@ -0,0 +1,221 @@
+function sparseimages(nm,typ)
+  if (strcmp(typ,"txt"))
+    txtimages(nm,15,typ);
+  else
+    if (strcmp (nm, "gplot"))
+      gplotimages("gplot",typ);
+    elseif (strcmp (nm, "grid"))
+      femimages("grid",typ);
+    else
+      otherimages(nm,200,typ);
+    endif
+  endif
+  ## Kluge to give gnuplot enough time to process last figure before we
+  ## exit.  Otherwise, Octave will delete the temporary data files when
+  ## it exits and gnuplot will fail...
+  sleep (1);
+endfunction
+
+## Use this function before plotting commands and after every call to
+## print since print() resets output to stdout (unfortunately, gnpulot
+## can't pop output as it can the terminal type).
+function bury_output ()
+  automatic_replot(0);
+  __gnuplot_set__ term dumb
+  [status, dummy] = fileattrib("/dev/null");
+  if (status)
+    __gnuplot_raw__ ("set output \"/dev/null\"\n");
+  endif
+endfunction
+
+function gplotimages(nm,typ)
+  bury_output ();
+  A = sparse([2,6,1,3,2,4,3,5,4,6,1,5],
+	     [1,1,2,2,3,3,4,4,5,5,6,6],1,6,6);
+  xy = [0,4,8,6,4,2;5,0,5,7,5,7]';
+  gplot(A,xy)
+  print(strcat(nm,".",typ),strcat("-d",typ))
+  bury_output ();
+endfunction
+
+function txtimages(nm,n,typ)
+  a = 10*speye(n) + sparse(1:n,ceil([1:n]/2),1,n,n) + ...
+      sparse(ceil([1:n]/2),1:n,1,n,n);
+  if (strcmp (nm, "gplot") || strcmp (nm, "grid"))
+    fid = fopen (sprintf ("%s.txt", nm), "wt");
+    fputs (fid, "+---------------------------------+\n");
+    fputs (fid, "| Image unavailable in text mode. |\n");
+    fputs (fid, "+---------------------------------+\n");
+    fclose (fid);
+  elseif (strcmp (nm, "spmatrix"))
+    printsparse(a,strcat("spmatrix.",typ));
+  else
+    if (!isempty(findstr(octave_config_info ("DEFS"),"HAVE_COLAMD")) &&
+	!isempty(findstr(octave_config_info ("DEFS"),"HAVE_CHOLMOD")))
+      if (strcmp (nm, "spchol"))
+	r1 = chol(a);
+	printsparse(r1,strcat("spchol.",typ));
+      elseif (strcmp (nm, "spcholperm"))
+	[r2,p2,q2]=chol(a);
+	printsparse(r2,strcat("spcholperm.",typ));
+      endif
+      ## printf("Text NNZ: Matrix %d, Chol %d, PermChol %d\n",nnz(a),nnz(r1),nnz(r2));
+    endif
+  endif
+endfunction
+
+function otherimages(nm,n,typ)
+  bury_output ();
+  a = 10*speye(n) + sparse(1:n,ceil([1:n]/2),1,n,n) + ...
+      sparse(ceil([1:n]/2),1:n,1,n,n);
+  if (strcmp (nm, "spmatrix"))
+    spy(a);
+    axis("ij")
+    print(strcat("spmatrix.",typ),strcat("-d",typ))
+    bury_output ();
+  else
+    if (!isempty(findstr(octave_config_info ("DEFS"),"HAVE_COLAMD")) &&
+	!isempty(findstr(octave_config_info ("DEFS"),"HAVE_CHOLMOD")))
+      if (strcmp (nm, "spchol"))
+	r1 = chol(a);
+	spy(r1);
+	axis("ij")
+	print(strcat("spchol.",typ),strcat("-d",typ))
+	bury_output ();
+      elseif (strcmp (nm, "spcholperm"))
+	[r2,p2,q2]=chol(a);
+	spy(r2);
+	axis("ij")
+	print(strcat("spcholperm.",typ),strcat("-d",typ))
+	bury_output ();
+      endif
+      ## printf("Image NNZ: Matrix %d, Chol %d, PermChol %d\n",nnz(a),nnz(r1),nnz(r2));
+    endif
+  endif
+endfunction
+
+function printsparse(a,nm)
+  fid = fopen (nm,"wt");
+  for i = 1:size(a,1)
+    if (rem(i,5) == 0)
+      fprintf (fid,"         %2d - ", i);
+    else
+      fprintf (fid,"            | ");
+    endif
+    for j = 1:size(a,2)
+      if (a(i,j) == 0)
+	fprintf(fid,"  ")
+      else
+	fprintf(fid," *")
+      endif
+    endfor
+    fprintf(fid,"\n")
+  endfor
+  fprintf(fid,"            |-");
+  for j=1:size(a,2)
+    if (rem(j,5)==0)
+      fprintf(fid,"-|");
+    else
+      fprintf(fid,"--");
+    endif
+  endfor
+  fprintf(fid,"\n")
+  fprintf(fid,"              ");
+  for j=1:size(a,2)
+    if (rem(j,5)==0)
+      fprintf(fid,"%2d",j);
+    else
+      fprintf(fid,"  ");
+    endif
+  endfor
+  fclose(fid);
+endfunction
+
+function femimages (nm,typ)
+  bury_output ();
+  if (!isempty(findstr(octave_config_info ("DEFS"),"HAVE_COLAMD")) &&
+      !isempty(findstr(octave_config_info ("DEFS"),"HAVE_CHOLMOD")) &&
+      !isempty(findstr(octave_config_info ("DEFS"),"HAVE_UMFPACK")))
+    ## build a rectangle
+    node_y = [1;1.2;1.5;1.8;2]*ones(1,11);
+    node_x = ones(5,1)*[1,1.05,1.1,1.2,1.3,1.5,1.7,1.8,1.9,1.95,2];
+    nodes = [node_x(:), node_y(:)];
+
+    [h,w] = size(node_x);
+    elems = [];
+    for idx = 1:w-1
+      widx = (idx-1)*h;
+      elems = [elems; widx+[(1:h-1);(2:h);h+(1:h-1)]']; 
+      elems = [elems; widx+[(2:h);h+(2:h);h+(1:h-1)]']; 
+    endfor
+
+    E = size(elems,1);  #No. of elements
+    N = size(nodes,1);  #No. of elements
+    D = size(elems,2);  #dimentions+1
+
+    ## Plot FEM Geometry
+    elemx = elems(:,[1,2,3,1])';
+    xelems = reshape( nodes(elemx, 1), 4, E);
+    yelems = reshape( nodes(elemx, 2), 4, E);
+
+    ## Set element conductivity
+    conductivity = [1*ones(1,16),2*ones(1,48),1*ones(1,16)];
+
+    ## Dirichlet boundary conditions
+    D_nodes = [1:5, 51:55]; 
+    D_value = [10*ones(1,5), 20*ones(1,5)]; 
+  
+    ## Neumann boundary conditions
+    ## Note that N_value must be normalized by the boundary
+    ##   length and element conductivity
+    N_nodes = [];
+    N_value = [];
+
+    ## Calculate connectivity matrix
+    C = sparse((1:D*E), reshape(elems',D*E,1),1, D*E, N);
+
+    ## Calculate stiffness matrix
+    Siidx = floor([0:D*E-1]'/D)*D*ones(1,D) + ones(D*E,1)*(1:D) ;
+    Sjidx = [1:D*E]'*ones(1,D);
+    Sdata = zeros(D*E,D);
+    dfact = prod(2:(D-1));
+    for j = 1:E
+      a = inv([ ones(D,1), nodes( elems(j,:), : ) ]);
+      const = conductivity(j)*2/dfact/abs(det(a));
+      Sdata(D*(j-1)+(1:D),:)= const * a(2:D,:)'*a(2:D,:);
+    endfor
+
+    ## Element-wise system matrix
+    SE = sparse(Siidx,Sjidx,Sdata);
+    ## Global system matrix
+    S = C'* SE *C;
+
+    ## Set Dirichlet boundary
+    V = zeros(N,1);
+    V(D_nodes) = D_value;
+    idx = 1:N;
+    idx(D_nodes) = [];
+
+    ## Set Neumann boundary
+    Q = zeros(N,1);
+    Q(N_nodes) = N_value; # FIXME
+
+    V(idx) = S(idx,idx)\( Q(idx) - S(idx,D_nodes)*V(D_nodes) );
+
+    velems = reshape( V(elemx), 4, E);
+
+    sz = size(xelems,2);
+    ## FIXME How can I do this without a gnuplot specific commands? plot3 anyone?
+    unwind_protect
+      __gnuplot_set__  parametric;
+      __gnuplot_raw__ ("set nohidden3d;\n");
+      tmp = [([xelems; NaN*ones(1,sz)])(:), ([yelems; NaN*ones(1,sz)])(:), ([velems; NaN*ones(1,sz)])(:)];
+      __gnuplot_splot__(tmp);
+      __gnuplot_raw__ ("set view 80,10;\n")
+      print(strcat(nm,".",typ),strcat("-d",typ))
+      bury_output ();
+    unwind_protect_cleanup
+      __gnuplot_set__ noparametric; 
+    end_unwind_protect
+  endif
+endfunction