--- a/doc/interpreter/sparse.txi	Thu Jun 01 16:16:01 2006 +0000
+++ b/doc/interpreter/sparse.txi	Thu Jun 01 19:05:32 2006 +0000
@@ -427,8 +427,8 @@
   sprank, svds, spaugment)
 
 @item Iterative techniques:
-  @dfn{luinc}, (bicg, bicgstab, cholinc, cgs, gmres, lsqr, minres, 
-  pcg, pcr, qmr, symmlq)
+  @dfn{luinc}, @dfn{pcg}, @dfn{pcr}, (bicg, bicgstab, cholinc, cgs, 
+  gmres, lsqr, minres, qmr, symmlq)
 
 @item Miscellaneous:
   @dfn{spparms}, @dfn{symbfact}, @dfn{spstats}, 
@@ -744,7 +744,10 @@
 @node Iterative Techniques, Real Life Example, Sparse Linear Algebra, Sparse Matrices
 @section Iterative Techniques applied to sparse matrices
 
-WRITE ME OR DELETE ME IF THERE ARE NO ITERATIVE SOLVERS IN OCTAVE 3.0
+There are three functions currently to document here, these being
+@dfn{luinc}, @dfn{pcg} and @dfn{pcr}.
+
+WRITE ME.
 
 @node Real Life Example, Oct-Files, Iterative Techniques, Sparse Matrices
 @section Real Life Example of the use of Sparse Matrices
@@ -1478,9 +1481,13 @@
 @item minres
 @emph{Not implemented}
 @item pcg
-@emph{Not implemented}
+Solves the linear system of equations @code{@var{A} * @var{x} =
+@var{b}} by means of the  Preconditioned Conjugate Gradient iterative
+method.
 @item pcr
-@emph{Not implemented}
+Solves the linear system of equations @code{@var{A} * @var{x} =
+@var{b}} by means of the  Preconditioned Conjugate Residual iterative
+method.
 @item qmr
 @emph{Not implemented}
 @item symmlq
@@ -1548,6 +1555,10 @@
 		matrix S
 * nzmax::	Returns the amount of storage allocated to the sparse
 		matrix SM.
+* pcg::		Solves linear system of equations by means of the 
+		Preconditioned Conjugate Gradient iterative method.
+* pcr::		Solves linear system of equations by means of the 
+		Preconditioned Conjugate Residual iterative method.
 * spalloc::	Returns an empty sparse matrix of size R-by-C.
 * sparse::	SPARSE: create a sparse matrix
 * spatan2::	Compute atan (Y / X) for corresponding sparse matrix
@@ -1672,12 +1683,22 @@
 
 @DOCSTRING(nonzeros)
 
-@node nzmax, spalloc, nonzeros, Function Reference
+@node nzmax, pcg, nonzeros, Function Reference
 @subsubsection nzmax
 
 @DOCSTRING(nzmax)
 
-@node spalloc, sparse, nzmax, Function Reference
+@node pcg, pcr, nzmax, Function Reference
+@subsubsection pcg
+
+@DOCSTRING(pcg)
+
+@node pcr, spalloc, pcg, Function Reference
+@subsubsection pcr
+
+@DOCSTRING(pcr)
+
+@node spalloc, sparse, pcr, Function Reference
 @subsubsection spalloc
 
 @DOCSTRING(spalloc)

--- a/libcruft/ChangeLog	Thu Jun 01 16:16:01 2006 +0000
+++ b/libcruft/ChangeLog	Thu Jun 01 19:05:32 2006 +0000
@@ -1,3 +1,7 @@
+2006-06-01  David Bateman  <dbateman@free.fr>
+
+	* slatec-fn/dpchim.f, slatec-fn/dpchst.f: New files.
+
 2006-05-22  John W. Eaton  <jwe@octave.org>
 
 	* lapack/dlantr.f, lapack/zlantr.f: New files.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/slatec-fn/dpchim.f	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,285 @@
+*DECK DPCHIM
+      SUBROUTINE DPCHIM (N, X, F, D, INCFD, IERR)
+C***BEGIN PROLOGUE  DPCHIM
+C***PURPOSE  Set derivatives needed to determine a monotone piecewise
+C            cubic Hermite interpolant to given data.  Boundary values
+C            are provided which are compatible with monotonicity.  The
+C            interpolant will have an extremum at each point where mono-
+C            tonicity switches direction.  (See DPCHIC if user control
+C            is desired over boundary or switch conditions.)
+C***LIBRARY   SLATEC (PCHIP)
+C***CATEGORY  E1A
+C***TYPE      DOUBLE PRECISION (PCHIM-S, DPCHIM-D)
+C***KEYWORDS  CUBIC HERMITE INTERPOLATION, MONOTONE INTERPOLATION,
+C             PCHIP, PIECEWISE CUBIC INTERPOLATION
+C***AUTHOR  Fritsch, F. N., (LLNL)
+C             Lawrence Livermore National Laboratory
+C             P.O. Box 808  (L-316)
+C             Livermore, CA  94550
+C             FTS 532-4275, (510) 422-4275
+C***DESCRIPTION
+C
+C          DPCHIM:  Piecewise Cubic Hermite Interpolation to
+C                  Monotone data.
+C
+C     Sets derivatives needed to determine a monotone piecewise cubic
+C     Hermite interpolant to the data given in X and F.
+C
+C     Default boundary conditions are provided which are compatible
+C     with monotonicity.  (See DPCHIC if user control of boundary con-
+C     ditions is desired.)
+C
+C     If the data are only piecewise monotonic, the interpolant will
+C     have an extremum at each point where monotonicity switches direc-
+C     tion.  (See DPCHIC if user control is desired in such cases.)
+C
+C     To facilitate two-dimensional applications, includes an increment
+C     between successive values of the F- and D-arrays.
+C
+C     The resulting piecewise cubic Hermite function may be evaluated
+C     by DPCHFE or DPCHFD.
+C
+C ----------------------------------------------------------------------
+C
+C  Calling sequence:
+C
+C        PARAMETER  (INCFD = ...)
+C        INTEGER  N, IERR
+C        DOUBLE PRECISION  X(N), F(INCFD,N), D(INCFD,N)
+C
+C        CALL  DPCHIM (N, X, F, D, INCFD, IERR)
+C
+C   Parameters:
+C
+C     N -- (input) number of data points.  (Error return if N.LT.2 .)
+C           If N=2, simply does linear interpolation.
+C
+C     X -- (input) real*8 array of independent variable values.  The
+C           elements of X must be strictly increasing:
+C                X(I-1) .LT. X(I),  I = 2(1)N.
+C           (Error return if not.)
+C
+C     F -- (input) real*8 array of dependent variable values to be
+C           interpolated.  F(1+(I-1)*INCFD) is value corresponding to
+C           X(I).  DPCHIM is designed for monotonic data, but it will
+C           work for any F-array.  It will force extrema at points where
+C           monotonicity switches direction.  If some other treatment of
+C           switch points is desired, DPCHIC should be used instead.
+C                                     -----
+C     D -- (output) real*8 array of derivative values at the data
+C           points.  If the data are monotonic, these values will
+C           determine a monotone cubic Hermite function.
+C           The value corresponding to X(I) is stored in
+C                D(1+(I-1)*INCFD),  I=1(1)N.
+C           No other entries in D are changed.
+C
+C     INCFD -- (input) increment between successive values in F and D.
+C           This argument is provided primarily for 2-D applications.
+C           (Error return if  INCFD.LT.1 .)
+C
+C     IERR -- (output) error flag.
+C           Normal return:
+C              IERR = 0  (no errors).
+C           Warning error:
+C              IERR.GT.0  means that IERR switches in the direction
+C                 of monotonicity were detected.
+C           "Recoverable" errors:
+C              IERR = -1  if N.LT.2 .
+C              IERR = -2  if INCFD.LT.1 .
+C              IERR = -3  if the X-array is not strictly increasing.
+C             (The D-array has not been changed in any of these cases.)
+C               NOTE:  The above errors are checked in the order listed,
+C                   and following arguments have **NOT** been validated.
+C
+C***REFERENCES  1. F. N. Fritsch and J. Butland, A method for construc-
+C                 ting local monotone piecewise cubic interpolants, SIAM
+C                 Journal on Scientific and Statistical Computing 5, 2
+C                 (June 1984), pp. 300-304.
+C               2. F. N. Fritsch and R. E. Carlson, Monotone piecewise
+C                 cubic interpolation, SIAM Journal on Numerical Ana-
+C                 lysis 17, 2 (April 1980), pp. 238-246.
+C***ROUTINES CALLED  DPCHST, XERMSG
+C***REVISION HISTORY  (YYMMDD)
+C   811103  DATE WRITTEN
+C   820201  1. Introduced  DPCHST  to reduce possible over/under-
+C             flow problems.
+C           2. Rearranged derivative formula for same reason.
+C   820602  1. Modified end conditions to be continuous functions
+C             of data when monotonicity switches in next interval.
+C           2. Modified formulas so end conditions are less prone
+C             of over/underflow problems.
+C   820803  Minor cosmetic changes for release 1.
+C   870707  Corrected XERROR calls for d.p. name(s).
+C   870813  Updated Reference 1.
+C   890206  Corrected XERROR calls.
+C   890411  Added SAVE statements (Vers. 3.2).
+C   890531  Changed all specific intrinsics to generic.  (WRB)
+C   890703  Corrected category record.  (WRB)
+C   890831  Modified array declarations.  (WRB)
+C   891006  Cosmetic changes to prologue.  (WRB)
+C   891006  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
+C   920429  Revised format and order of references.  (WRB,FNF)
+C***END PROLOGUE  DPCHIM
+C  Programming notes:
+C
+C     1. The function  DPCHST(ARG1,ARG2)  is assumed to return zero if
+C        either argument is zero, +1 if they are of the same sign, and
+C        -1 if they are of opposite sign.
+C     2. To produce a single precision version, simply:
+C        a. Change DPCHIM to PCHIM wherever it occurs,
+C        b. Change DPCHST to PCHST wherever it occurs,
+C        c. Change all references to the Fortran intrinsics to their
+C           single precision equivalents,
+C        d. Change the double precision declarations to real, and
+C        e. Change the constants ZERO and THREE to single precision.
+C
+C  DECLARE ARGUMENTS.
+C
+      INTEGER  N, INCFD, IERR
+      DOUBLE PRECISION  X(*), F(INCFD,*), D(INCFD,*)
+C
+C  DECLARE LOCAL VARIABLES.
+C
+      INTEGER  I, NLESS1
+      DOUBLE PRECISION  DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE,
+     *      H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO
+      SAVE ZERO, THREE
+      DOUBLE PRECISION  DPCHST
+      DATA  ZERO /0.D0/, THREE/3.D0/
+C
+C  VALIDITY-CHECK ARGUMENTS.
+C
+C***FIRST EXECUTABLE STATEMENT  DPCHIM
+      IF ( N.LT.2 )  GO TO 5001
+      IF ( INCFD.LT.1 )  GO TO 5002
+      DO 1  I = 2, N
+         IF ( X(I).LE.X(I-1) )  GO TO 5003
+    1 CONTINUE
+C
+C  FUNCTION DEFINITION IS OK, GO ON.
+C
+      IERR = 0
+      NLESS1 = N - 1
+      H1 = X(2) - X(1)
+      DEL1 = (F(1,2) - F(1,1))/H1
+      DSAVE = DEL1
+C
+C  SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION.
+C
+      IF (NLESS1 .GT. 1)  GO TO 10
+      D(1,1) = DEL1
+      D(1,N) = DEL1
+      GO TO 5000
+C
+C  NORMAL CASE  (N .GE. 3).
+C
+   10 CONTINUE
+      H2 = X(3) - X(2)
+      DEL2 = (F(1,3) - F(1,2))/H2
+C
+C  SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
+C     SHAPE-PRESERVING.
+C
+      HSUM = H1 + H2
+      W1 = (H1 + HSUM)/HSUM
+      W2 = -H1/HSUM
+      D(1,1) = W1*DEL1 + W2*DEL2
+      IF ( DPCHST(D(1,1),DEL1) .LE. ZERO)  THEN
+         D(1,1) = ZERO
+      ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO)  THEN
+C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
+         DMAX = THREE*DEL1
+         IF (ABS(D(1,1)) .GT. ABS(DMAX))  D(1,1) = DMAX
+      ENDIF
+C
+C  LOOP THROUGH INTERIOR POINTS.
+C
+      DO 50  I = 2, NLESS1
+         IF (I .EQ. 2)  GO TO 40
+C
+         H1 = H2
+         H2 = X(I+1) - X(I)
+         HSUM = H1 + H2
+         DEL1 = DEL2
+         DEL2 = (F(1,I+1) - F(1,I))/H2
+   40    CONTINUE
+C
+C        SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC.
+C
+         D(1,I) = ZERO
+         IF ( DPCHST(DEL1,DEL2) .LT. 0.)  GO TO 42
+         IF ( DPCHST(DEL1,DEL2) .EQ. 0.)  GO TO 41
+         GO TO 45
+C
+C        COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY.
+C
+   41    CONTINUE
+         IF (DEL2 .EQ. ZERO)  GO TO 50
+         IF ( DPCHST(DSAVE,DEL2) .LT. ZERO)  IERR = IERR + 1
+         DSAVE = DEL2
+         GO TO 50
+C
+   42    CONTINUE
+         IERR = IERR + 1
+         DSAVE = DEL2
+         GO TO 50
+C
+C        USE BRODLIE MODIFICATION OF BUTLAND FORMULA.
+C
+   45    CONTINUE
+         HSUMT3 = HSUM+HSUM+HSUM
+         W1 = (HSUM + H1)/HSUMT3
+         W2 = (HSUM + H2)/HSUMT3
+         DMAX = MAX( ABS(DEL1), ABS(DEL2) )
+         DMIN = MIN( ABS(DEL1), ABS(DEL2) )
+         DRAT1 = DEL1/DMAX
+         DRAT2 = DEL2/DMAX
+         D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2)
+C
+   50 CONTINUE
+C
+C  SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
+C     SHAPE-PRESERVING.
+C
+      W1 = -H2/HSUM
+      W2 = (H2 + HSUM)/HSUM
+      D(1,N) = W1*DEL1 + W2*DEL2
+      IF ( DPCHST(D(1,N),DEL2) .LE. ZERO)  THEN
+         D(1,N) = ZERO
+      ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO)  THEN
+C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
+         DMAX = THREE*DEL2
+         IF (ABS(D(1,N)) .GT. ABS(DMAX))  D(1,N) = DMAX
+      ENDIF
+C
+C  NORMAL RETURN.
+C
+ 5000 CONTINUE
+      RETURN
+C
+C  ERROR RETURNS.
+C
+ 5001 CONTINUE
+C     N.LT.2 RETURN.
+      IERR = -1
+      CALL XERMSG ('SLATEC', 'DPCHIM',
+     +   'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
+      RETURN
+C
+ 5002 CONTINUE
+C     INCFD.LT.1 RETURN.
+      IERR = -2
+      CALL XERMSG ('SLATEC', 'DPCHIM', 'INCREMENT LESS THAN ONE', IERR,
+     +   1)
+      RETURN
+C
+ 5003 CONTINUE
+C     X-ARRAY NOT STRICTLY INCREASING.
+      IERR = -3
+      CALL XERMSG ('SLATEC', 'DPCHIM',
+     +   'X-ARRAY NOT STRICTLY INCREASING', IERR, 1)
+      RETURN
+C------------- LAST LINE OF DPCHIM FOLLOWS -----------------------------
+      END

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/slatec-fn/dpchst.f	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,59 @@
+*DECK DPCHST
+      DOUBLE PRECISION FUNCTION DPCHST (ARG1, ARG2)
+C***BEGIN PROLOGUE  DPCHST
+C***SUBSIDIARY
+C***PURPOSE  DPCHIP Sign-Testing Routine
+C***LIBRARY   SLATEC (PCHIP)
+C***TYPE      DOUBLE PRECISION (PCHST-S, DPCHST-D)
+C***AUTHOR  Fritsch, F. N., (LLNL)
+C***DESCRIPTION
+C
+C         DPCHST:  DPCHIP Sign-Testing Routine.
+C
+C
+C     Returns:
+C        -1. if ARG1 and ARG2 are of opposite sign.
+C         0. if either argument is zero.
+C        +1. if ARG1 and ARG2 are of the same sign.
+C
+C     The object is to do this without multiplying ARG1*ARG2, to avoid
+C     possible over/underflow problems.
+C
+C  Fortran intrinsics used:  SIGN.
+C
+C***SEE ALSO  DPCHCE, DPCHCI, DPCHCS, DPCHIM
+C***ROUTINES CALLED  (NONE)
+C***REVISION HISTORY  (YYMMDD)
+C   811103  DATE WRITTEN
+C   820805  Converted to SLATEC library version.
+C   870813  Minor cosmetic changes.
+C   890411  Added SAVE statements (Vers. 3.2).
+C   890531  Changed all specific intrinsics to generic.  (WRB)
+C   890531  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900328  Added TYPE section.  (WRB)
+C   910408  Updated AUTHOR and DATE WRITTEN sections in prologue.  (WRB)
+C   930503  Improved purpose.  (FNF)
+C***END PROLOGUE  DPCHST
+C
+C**End
+C
+C  DECLARE ARGUMENTS.
+C
+      DOUBLE PRECISION  ARG1, ARG2
+C
+C  DECLARE LOCAL VARIABLES.
+C
+      DOUBLE PRECISION  ONE, ZERO
+      SAVE ZERO, ONE
+      DATA  ZERO /0.D0/,  ONE/1.D0/
+C
+C  PERFORM THE TEST.
+C
+C***FIRST EXECUTABLE STATEMENT  DPCHST
+      DPCHST = SIGN(ONE,ARG1) * SIGN(ONE,ARG2)
+      IF ((ARG1.EQ.ZERO) .OR. (ARG2.EQ.ZERO))  DPCHST = ZERO
+C
+      RETURN
+C------------- LAST LINE OF DPCHST FOLLOWS -----------------------------
+      END

--- a/liboctave/Array.cc	Thu Jun 01 16:16:01 2006 +0000
+++ b/liboctave/Array.cc	Thu Jun 01 19:05:32 2006 +0000
@@ -3104,14 +3104,15 @@
 	      else
 		final_lhs_dims = new_dims;
 
-	      lhs_dims = new_dims;
-
-	      lhs_dims_len = lhs_dims.length ();
-
-	      frozen_len = freeze (idx, lhs_dims, true);
+	      lhs_dims_len = new_dims.length ();
+
+	      frozen_len = freeze (idx, new_dims, true);
 
 	      if (rhs_is_scalar)
 		{
+		  if (n_idx < orig_lhs_dims_len)
+		    lhs = lhs.reshape (lhs_dims);
+
 		  lhs.resize_and_fill (new_dims, rfv);
 
 		  if  (! final_lhs_dims.any_zero ())
@@ -3147,6 +3148,9 @@
 		    }
 		  else
 		    {
+		      if (n_idx < orig_lhs_dims_len)
+			lhs = lhs.reshape (lhs_dims);
+
 		      lhs.resize_and_fill (new_dims, rfv);
 
 		      if  (! final_lhs_dims.any_zero ())

--- a/liboctave/ChangeLog	Thu Jun 01 16:16:01 2006 +0000
+++ b/liboctave/ChangeLog	Thu Jun 01 19:05:32 2006 +0000
@@ -1,3 +1,7 @@
+2006-05-31  David Bateman  <dbateman@free.fr>
+
+	* Array.cc (assignN): Maybe reshape LHS before doing assignment.
+
 2006-05-23  John W. Eaton  <jwe@octave.org>
 
 	* oct-types.h.in: Include stdint.h or inttypes.h for integer

--- a/scripts/ChangeLog	Thu Jun 01 16:16:01 2006 +0000
+++ b/scripts/ChangeLog	Thu Jun 01 19:05:32 2006 +0000
@@ -1,3 +1,11 @@
+2006-06-01  David Bateman  <dbateman@free.fr>
+
+	* general/interpft.m, general/quadl.m, general/polyarea.m,
+	general/interp1.m, general/gradient.m, general/interp2.m,
+	general/bicubic.m, miscellaneous/inputname.m, plot/__plt3__.m,
+	plot/ndgrid.m, plot/plot3.m, polynomial/pchip.m, sparse/pcg.m,
+	sparse/pcr.m, strings/mat2str.m: New files from Octave Forge.
+
 2006-05-31  Bill Denney  <bill@givebillmoney.com>
 
  	* miscellaneous/fileparts.m, miscellaneous/fullfile.m: Add seealso.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/bicubic.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,197 @@
+## Copyright (C) 2005  Hoxide Ma
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{zi}=} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
+##
+## Bicubic interpolation method. Returns a matrix @var{zi}, corresponding the
+## the bicubic interpolations at @var{xi} and @var{yi} of the data supplied
+## as @var{x}, @var{y} and @var{z}. 
+##
+## For further information please see bicubic.pdf available at
+## @url{http://wiki.woodpecker.org.cn/moin/Octave/Bicubic}
+## @end deftypefn
+## @seealso{interp2}
+
+## Bicubic interpolation method.
+## Author: Hoxide Ma <hoxide_dirac@yahoo.com.cn>
+
+function F = bicubic(X, Y, Z, XI, YI, spline_alpha)
+
+  if (nargin < 1 || nargin > 6)
+    print_usage ();
+  endif
+
+  if (nargin == 6 && prod(size(spline_alpha))==1)
+    a = spline_alpha
+  else
+    a = 0.5;
+  endif
+
+  if (nargin <=2) # bicubic(Z) or bicubic(Z, 2)
+    if (nargin == 1) 
+      n = 1;
+    else
+      n = Y;
+    endif
+    Z= X;
+    X = [];
+    [rz,cz] = size(Z);
+    s = linspace(1, cz, (cz-1)*pow2(n)+1);
+    t = linspace(1, rz, (rz-1)*pow2(n)+1);
+  elseif (nargin == 3)
+    if (!isvector (X) || !isvector (Y))
+      error ("XI and YI must be vector");
+    endif
+    s = Y;
+    t = Z;
+    Z = X;
+    [rz, cz] = size(Z);
+  elseif (nargin == 5 || nargin == 6)
+    [rz, cz] = size (Z) ; 
+    if (isvector (X) && isvector (Y))
+      if (rz != length (Y) || cz != length (X))
+	error ("length of X and Y must match the size of Z");
+      endif
+    elseif (size(X) == size(Y) && size(X) == size(Z))
+      X = X(1,:);
+      Y = Y(:,1);
+    else
+      error ("X, Y and Z must be martrices of same size");
+    endif
+    
+    ## mark values outside the lookup table
+    xfirst_ind = find (XI < X(1));
+    xlast_ind  = find (XI > X(cz));    
+    yfirst_ind = find (YI < Y(1));
+    ylast_ind  = find (YI > Y(rz));
+    ## set value outside the table preliminary to min max index   
+    XI(xfirst_ind) = X(1);
+    XI(xlast_ind) = X(cz);
+    YI(yfirst_ind) = Y(1);
+    YI(ylast_ind) = Y(rz);
+
+
+    X = reshape(X, 1, cz);
+    X(cz) *= (1+(sign(X(cz)))*eps);
+    if (X(cz) == 0) 
+      X(cz) = eps;
+    endif; 
+    XI = reshape(XI,1,length(XI));
+    [m, i] = sort([X, XI]);
+    o = cumsum ( i<= cz);
+    xidx = o([find( i> cz)]);
+    
+    Y = reshape(Y, rz, 1);
+    Y(rz) *= (1+(sign(Y(rz)))*eps);
+    if (Y(rz) == 0) 
+      Y(rz) = eps;
+    endif; 
+    YI = reshape(YI,length(YI),1);
+    [m, i] = sort([Y; YI]);
+    o = cumsum ( i<= rz);
+    yidx = o([find( i> rz)]);
+    
+    ## set s and t used follow codes
+    s = xidx + ((XI .- X(xidx))./(X(xidx+1) .- X(xidx)));
+    t = yidx + ((YI - Y(yidx))./(Y(yidx+1) - Y(yidx)));
+  else
+    print_usage ();
+  endif
+  
+  if (rz < 3 || cz < 3)
+    error ("Z at least a 3 by 3 matrices");
+  endif
+
+  inds = floor(s);
+  d = find(s==cz);
+  s = s - floor(s);
+  inds(d) = cz-1;
+  s(d) = 1.0;
+  
+  d = [];
+  indt = floor(t);
+  d = find(t==rz);
+  t = t - floor(t);
+  indt(d) = rz-1;
+  t(d) = 1.0;
+  d = [];
+
+  p = zeros(size(Z)+2);
+  p(2:rz+1,2:cz+1) = Z;
+  p(1,:)      =    (6*(1-a))*p(2,:)    -3*p(3,:)   + (6*a-2)*p(4,:);
+  p(rz+2,:)   =    (6*(1-a))*p(rz+1,:) -3*p(rz,:)  + (6*a-2)*p(rz-1,:);
+  p(:,1)      =    (6*(1-a))*p(:,2)    -3*p(:,3)   + (6*a-2)*p(:,4);
+  p(:,cz+2)   =    (6*(1-a))*p(:,cz+1)  -3*p(:,cz) + (6*a-2)*p(:,cz-1);
+
+  ## calculte the C1(t) C2(t) C3(t) C4(t) and C1(s) C2(s) C3(s) C4(s)
+  t2= t.*t;
+  t3= t2.*t;
+
+  ct0=     -a .* t3 +     (2 * a) .* t2 - a .* t ;      # -a G0
+  ct1 = (2-a) .* t3 +      (-3+a) .* t2          + 1 ;  # F0 - a G1
+  ct2 = (a-2) .* t3 + (-2 *a + 3) .* t2 + a .* t ;      # F1 + a G0
+  ct3 =     a .* t3 -           a .* t2;                # a G1
+  t = [];t2=[]; t3=[];
+
+  s2= s.*s;
+  s3= s2.*s;
+
+  cs0=     -a .* s3 +     (2 * a) .* s2 - a .*s ;      # -a G0
+  cs1 = (2-a) .* s3 +    (-3 + a) .* s2         + 1 ;  # F0 - a G1
+  cs2 = (a-2) .* s3 + (-2 *a + 3) .* s2 + a .*s ;      # F1 + a G0
+  cs3 =     a .* s3 -           a .* s2;               # a G1
+  s=[] ; s2 = []; s3 = [];
+
+  cs0 = cs0([1,1,1,1],:);
+  cs1 = cs1([1,1,1,1],:);
+  cs2 = cs2([1,1,1,1],:);
+  cs3 = cs3([1,1,1,1],:);
+
+  lent = length(ct0);
+  lens = length(cs0);
+  F = zeros(lent,lens);
+  
+  for i = 1:lent
+    it = indt(i);
+    int = [it, it+1, it+2, it+3];
+    F(i,:) = [ct0(i),ct1(i),ct2(i),ct3(i)] * ...
+	(p(int,inds) .* cs0 + p(int, inds+1) .* cs1 + ...
+	 p(int, inds+2) .* cs2 + p(int, inds+3) .* cs3);
+  endfor
+
+  ## set points outside the table to NaN
+  if (! (isempty (xfirst_ind) && isempty (xlast_ind)))
+    F(:, [xfirst_ind, xlast_ind]) = NaN;
+  endif
+  if (! (isempty (yfirst_ind) && isempty (ylast_ind)))
+    F([yfirst_ind; ylast_ind], :) = NaN;
+  endif
+
+endfunction
+
+%!demo
+%! A=[13,-1,12;5,4,3;1,6,2];
+%! x=[0,1,4]+10; y=[-10,-9,-8];
+%! xi=linspace(min(x),max(x),17);
+%! yi=linspace(min(y),max(y),26);
+%! mesh(xi,yi,bicubic(x,y,A,xi,yi));
+%! [x,y] = meshgrid(x,y);
+%! __gnuplot_raw__ ("set nohidden3d;\n")
+%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/gradient.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,132 @@
+## Copyright (C) 2000  Kai Habel
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{x} = } gradient (@var{M})
+## @deftypefnx {Function File} {[@var{x}, @var{y}, @dots{}] = } gradient (@var{M})
+## @deftypefnx {Function File} {[@dots{}] = } gradient (@var{M}, @var{s})
+## @deftypefnx {Function File} {[@dots{}] = } gradient (@var{M}, @var{dx}, @var{dy}, @dots{})
+##
+## Calculates the gradient. @code{@var{x} = gradient (@var{M})}
+## calculates the one dimensional gradient if @var{M} is a vector. If
+## @var{M} is a matrix the gradient is calculated for each row.
+##
+## @code{[@var{x}, @var{y}] = gradient (@var{M})} calculates the one
+## dimensional gradient for each direction if @var{M} if @var{M} is a
+## matrix. Additional return arguments can be use for multi-dimensional
+## matrices.
+##
+## Spacing values between two points can be provided by the
+## @var{dx}, @var{dy} or @var{h} parameters. If @var{h} is supplied it
+## is assumed to be the spacing in all directions. Otherwise, seperate
+## values of the spacing can be supplied by the @var{dx}, etc variables.
+## A scalar value specifies an equidistant spacing, while a vector value
+## can be used to specify a variable spacing. The length must match
+## their respective dimension of @var{M}.
+## 
+## At boundary points a linear extrapolation is applied. Interior points
+## are calculated with the first approximation of the numerical gradient
+##
+## @example
+## y'(i) = 1/(x(i+1)-x(i-1)) *(y(i-1)-y(i+1)).
+## @end example
+## 
+## @end deftypefn
+
+## Author:  Kai Habel <kai.habel@gmx.de>
+## Modified: David Bateman <dbateman@free.fr> Added NDArray support
+
+function [varargout] = gradient (M, varargin)
+  
+  if ((nargin < 1))
+    print_usage ()
+  endif
+
+  transposed = false;
+  if (isvector (M))
+    ## make a column vector
+    transposed = (size(M,2) == 1);
+    M = M(:)';
+  endif
+
+  nd = ndims (M);
+  sz = size (M);
+  if (nargin > 2 && nargin != nd + 1)
+    print_usage ()
+  endif
+    
+  d = cell(1,nd);
+  if (nargin == 1)
+    for i=1:nd
+      d{i} = ones(sz(i), 1);
+    endfor
+  elseif (nargin == 2)
+    if (isscalar (varargin{1}))
+      for i=1:nd
+	d{i} = varargin{1} * ones(sz(i), 1);
+      endfor
+    else
+      for i=1:nd
+	d{i} = varargin{1};
+      endfor
+    endif
+  else
+    for i=1:nd
+      if (isscalar (varargin{1}))
+	d{i} = varargin{i} * ones(sz(i), 1);
+      else
+	d{i} = varargin{i};
+      endif
+    endfor
+
+    ## Why the hell did matlab decide to swap these two values?
+    tmp = d{1};
+    d{1} = d{2};
+    d{2} = tmp;
+  endif
+
+  for i = 1:max(2,min(nd,nargout))
+    mr = sz(i);
+    mc = prod([sz(1:i-1),sz(i+1:nd)]);
+    Y = zeros (size(M), class(M));
+
+    if (mr > 1)
+      ## top and bottom boundary
+      Y(1, :) = diff (M(1:2, :)) / d{i}(1);
+      Y(mr, :) = diff (M(mr - 1:mr, :)) / d{i}(mr - 1);
+    endif
+
+    if (mr > 2)
+      ## interior points
+      Y(2:mr-1, :) = (M(3:mr, :) .- M(1:mr - 2, :)) ./ ...
+          kron (d{i}(1:mr - 2) .+ d{i}(2:mr - 1), ones(1, mc));
+    endif
+    varargout{i} = ipermute (Y, [i:nd,1:i-1]);
+    M = permute (M, [2:nd,1]);
+  endfor
+
+  ## Why the hell did matlab decide to swap these two values?
+  tmp = varargout{1};
+  varargout{1} = varargout{2};
+  varargout{2} = tmp;
+
+  if (transposed)
+    varargout{1} = varargout{1}.';
+  endif
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/interp1.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,442 @@
+## Copyright (C) 2000 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{yi} =} interp1 (@var{x}, @var{y}, @var{xi})
+## @deftypefnx {Function File} {@var{yi} =} interp1 (@dots{}, @var{method})
+## @deftypefnx {Function File} {@var{yi} =} interp1 (@dots{}, @var{extrap})
+## @deftypefnx {Function File} {@var{pp} =} interp1 (@dots{}, 'pp')
+##
+## One-dimensional interpolation. Interpolate @var{y}, defined at the
+## points @var{x}, at the points @var{xi}. The sample points @var{x} 
+## must be strictly monotonic. If @var{y} is an array, treat the columns
+## of @var{y} seperately.
+##
+## Method is one of:
+##
+## @table @asis
+## @item 'nearest'
+## Return the nearest neighbour.
+## @item 'linear'
+## Linear interpolation from nearest neighbours
+## @item 'pchip'
+## Piece-wise cubic hermite interpolating polynomial
+## @item 'cubic'
+## Cubic interpolation from four nearest neighbours
+## @item 'spline'
+## Cubic spline interpolation--smooth first and second derivatives
+## throughout the curve
+## @end table
+##
+## Appending '*' to the start of the above method forces @code{interp1}
+## to assume that @var{x} is uniformly spaced, and only @code{@var{x}
+## (1)} and @code{@var{x} (2)} are referenced. This is usually faster,
+## and is never slower. The default method is 'linear'.
+##
+## If @var{extrap} is the string 'extrap', then extrapolate values beyond
+## the endpoints.  If @var{extrap} is a number, replace values beyond the
+## endpoints with that number.  If @var{extrap} is missing, assume NaN.
+##
+## If the string argument 'pp' is specified, then @var{xi} should not be
+## supplied and @code{interp1} returns the piece-wise polynomial that
+## can later be used with @code{ppval} to evaluate the interpolation.
+## There is an equivalence, such that @code{ppval (interp1 (@var{x},
+## @var{y}, @var{method}, 'pp'), @var{xi}) == interp1 (@var{x}, @var{y},
+## @var{xi}, @var{method}, 'extrap')}.
+##
+## An example of the use of @code{interp1} is
+##
+## @example
+## @group
+##    xf=[0:0.05:10]; yf = sin(2*pi*xf/5);
+##    xp=[0:10];      yp = sin(2*pi*xp/5);
+##    lin=interp1(xp,yp,xf);
+##    spl=interp1(xp,yp,xf,'spline');
+##    cub=interp1(xp,yp,xf,'cubic');
+##    near=interp1(xp,yp,xf,'nearest');
+##    plot(xf,yf,';original;',xf,lin,';linear;',xf,spl,';spline;',...
+##         xf,cub,';cubic;',xf,near,';nearest;',xp,yp,'*;;');
+## @end group
+## @end example
+##
+## @seealso{interpft}
+## @end deftypefn
+
+## 2000-03-25 Paul Kienzle
+##    added 'nearest' as suggested by Kai Habel
+## 2000-07-17 Paul Kienzle
+##    added '*' methods and matrix y
+##    check for proper table lengths
+## 2002-01-23 Paul Kienzle
+##    fixed extrapolation
+
+function yi = interp1(x, y, varargin)
+
+  if ( nargin < 3 || nargin > 6)
+    print_usage ();
+  endif
+
+  method = "linear";
+  extrap = NaN;
+  xi = [];
+  pp = false;
+  firstnumeric = true;
+
+  if (nargin > 2)
+    for i = 1:length(varargin)
+      arg = varargin{i};
+      if (ischar(arg))
+	arg = tolower (arg);
+	if (strcmp("extrap",arg))
+	  extrap = "extrap";
+	elseif (strcmp("pp",arg))
+	  pp = true;
+	else
+	  method = arg;
+	endif
+      else
+	if (firstnumeric)
+	  xi = arg;
+	  firstnumeric = false;
+	else
+	  extrap = arg;
+	endif
+      endif
+    endfor
+  endif
+
+  ## reshape matrices for convenience
+  x = x(:);
+  nx = size(x,1);
+  if (isvector(y) && size (y, 1) == 1)
+    y = y (:);
+  endif
+  ndy = ndims (y);
+  szy = size(y);
+  ny = szy(1);
+  nc = prod (szy(2:end));
+  y = reshape (y, ny, nc);
+  szx = size(xi);
+  xi = xi(:);
+
+  ## determine sizes
+  if (nx < 2 || ny < 2)
+     error ("interp1: table too short");
+  endif
+
+  ## determine which values are out of range and set them to extrap,
+  ## unless extrap=="extrap" in which case, extrapolate them like we
+  ## should be doing in the first place.
+  minx = x(1);
+  if (method(1) == "*")
+     dx = x(2) - x(1);
+     maxx = minx + (ny-1)*dx;
+  else
+     maxx = x(nx);
+  endif
+  if (!pp)
+    if ischar(extrap) && strcmp(extrap,"extrap")
+      range=1:size(xi,1);
+      yi = zeros(size(xi,1), size(y,2));
+    else
+      range = find(xi >= minx & xi <= maxx);
+      yi = extrap*ones(size(xi,1), size(y,2));
+      if isempty(range), 
+	if (!isvector(y) && length(szx) == 2 && (szx(1) == 1 || szx(2) == 1))
+	  if (szx(1) == 1)
+	    yi = reshape (yi, [szx(2), szy(2:end)]);
+	  else
+	    yi = reshape (yi, [szx(1), szy(2:end)]);
+	  endif
+	else
+	  yi = reshape (yi, [szx, szy(2:end)]);
+        endif
+        return; 
+      endif
+      xi = xi(range);
+    endif
+  endif
+
+  if strcmp(method, "nearest")
+    if (pp)
+      yi = mkpp ([x(1);(x(1:end-1)+x(2:end))/2;x(end)],y,szy(2:end));
+    else
+      idx = lookup(0.5*(x(1:nx-1)+x(2:nx)), xi)+1;
+      yi(range,:) = y(idx,:);
+    endif
+  elseif strcmp(method, "*nearest")
+    if (pp)
+      yi = mkpp ([minx; minx + [0.5:(ny-1)]'*dx; maxx],y,szy(2:end));
+    else
+      idx = max(1,min(ny,floor((xi-minx)/dx+1.5)));
+      yi(range,:) = y(idx,:);
+    endif
+  elseif strcmp(method, "linear")
+    dy = y(2:ny,:) - y(1:ny-1,:);
+    dx = x(2:nx) - x(1:nx-1);
+    if (pp)
+      yi = mkpp(x, [dy./dx, y(1:end-1)],szy(2:end));
+    else
+      ## find the interval containing the test point
+      idx = lookup (x(2:nx-1), xi)+1; 
+				# 2:n-1 so that anything beyond the ends
+				# gets dumped into an interval
+      ## use the endpoints of the interval to define a line
+      s = (xi - x(idx))./dx(idx);
+      yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:);
+    endif
+  elseif strcmp(method, "*linear")
+    if (pp)
+      dy = [y(2:ny,:) - y(1:ny-1,:)];
+      yi = mkpp (minx + [0:ny-1]*dx, [dy ./ dx, y(1:end-1)], szy(2:end));
+    else
+      ## find the interval containing the test point
+      t = (xi - minx)/dx + 1;
+      idx = max(1,min(ny,floor(t)));
+
+      ## use the endpoints of the interval to define a line
+      dy = [y(2:ny,:) - y(1:ny-1,:); y(ny,:) - y(ny-1,:)];
+      s = t - idx;
+      yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:); 
+    endif
+  elseif strcmp(method, "pchip") || strcmp(method, "*pchip")
+    if (nx == 2 || method(1) == "*") 
+      x = linspace(minx, maxx, ny); 
+    endif
+    ## Note that pchip's arguments are transposed relative to interp1
+    if (pp)
+      yi = pchip(x.', y.');
+      yi.d = szy(2:end);
+    else
+      yi(range,:) = pchip(x.', y.', xi.').';
+    endif
+
+  elseif strcmp(method, "cubic") || (strcmp(method, "*cubic") && pp)
+    ## FIXME Is there a better way to treat pp return return and *cubic
+    if (method(1) == "*") 
+      x = linspace(minx, maxx, ny).'; 
+      nx = ny;
+    endif
+
+    if (nx < 4 || ny < 4)
+      error ("interp1: table too short");
+    endif
+    idx = lookup(x(3:nx-2), xi) + 1;
+
+    ## Construct cubic equations for each interval using divided
+    ## differences (computation of c and d don't use divided differences
+    ## but instead solve 2 equations for 2 unknowns). Perhaps
+    ## reformulating this as a lagrange polynomial would be more efficient.
+    i=1:nx-3;
+    J = ones(1,nc);
+    dx = diff(x);
+    dx2 = x(i+1).^2 - x(i).^2;
+    dx3 = x(i+1).^3 - x(i).^3;
+    a=diff(y,3)./dx(i,J).^3/6;
+    b=(diff(y(1:nx-1,:),2)./dx(i,J).^2 - 6*a.*x(i+1,J))/2;
+    c=(diff(y(1:nx-2,:),1) - a.*dx3(:,J) - b.*dx2(:,J))./dx(i,J);
+    d=y(i,:) - ((a.*x(i,J) + b).*x(i,J) + c).*x(i,J);
+
+    if (pp)
+      xs = [x(1);x(3:nx-2)];
+      yi = mkpp ([x(1);x(3:nx-2);x(nx)], 
+		 [a(:), (b(:) + 3.*xs(:,J).*a(:)), ... 
+		  (c(:) + 2.*xs(:,J).*b(:) + 3.*xs(:,J)(:).^2.*a(:)), ...
+		  (d(:) + xs(:,J).*c(:) + xs(:,J).^2.*b(:) + ...
+		   xs(:,J).^3.*a(:))], szy(2:end));
+    else
+      yi(range,:) = ((a(idx,:).*xi(:,J) + b(idx,:)).*xi(:,J) ...
+		     + c(idx,:)).*xi(:,J) + d(idx,:);
+    endif
+  elseif strcmp(method, "*cubic")
+    if (nx < 4 || ny < 4)
+      error ("interp1: table too short");
+    endif
+
+    ## From: Miloje Makivic 
+    ## http://www.npac.syr.edu/projects/nasa/MILOJE/final/node36.html
+    t = (xi - minx)/dx + 1;
+    idx = max(min(floor(t), ny-2), 2);
+    t = t - idx;
+    t2 = t.*t;
+    tp = 1 - 0.5*t;
+    a = (1 - t2).*tp;
+    b = (t2 + t).*tp;
+    c = (t2 - t).*tp/3;
+    d = (t2 - 1).*t/6;
+    J = ones(1,nc);
+
+    yi(range,:) = a(:,J) .* y(idx,:) + b(:,J) .* y(idx+1,:) ...
+		  + c(:,J) .* y(idx-1,:) + d(:,J) .* y(idx+2,:);
+
+  elseif strcmp(method, "spline") || strcmp(method, "*spline")
+    if (nx == 2 || method(1) == "*") 
+      x = linspace(minx, maxx, ny); 
+    endif
+    ## Note that spline's arguments are transposed relative to interp1
+    if (pp)
+      yi = spline(x.', y.');
+      yi.d = szy(2:end);
+    else
+      yi(range,:) = spline(x.', y.', xi.').';
+    endif
+  else
+    error(["interp1 doesn't understand method '", method, "'"]);
+  endif
+
+  if (!pp)
+    if (!isvector(y) && length(szx) == 2 && (szx(1) == 1 || szx(2) == 1))
+      if (szx(1) == 1)
+	yi = reshape (yi, [szx(2), szy(2:end)]);
+      else
+	yi = reshape (yi, [szx(1), szy(2:end)]);
+      endif
+    else
+      yi = reshape (yi, [szx, szy(2:end)]);
+    endif
+  endif
+
+endfunction
+
+%!demo
+%! xf=0:0.05:10; yf = sin(2*pi*xf/5);
+%! xp=0:10;      yp = sin(2*pi*xp/5);
+%! lin=interp1(xp,yp,xf,"linear");
+%! spl=interp1(xp,yp,xf,"spline");
+%! cub=interp1(xp,yp,xf,"pchip");
+%! near=interp1(xp,yp,xf,"nearest");
+%! plot(xf,yf,";original;",xf,near,";nearest;",xf,lin,";linear;",...
+%!      xf,cub,";pchip;",xf,spl,";spline;",xp,yp,"*;;");
+%! %--------------------------------------------------------
+%! % confirm that interpolated function matches the original
+
+%!demo
+%! xf=0:0.05:10; yf = sin(2*pi*xf/5);
+%! xp=0:10;      yp = sin(2*pi*xp/5);
+%! lin=interp1(xp,yp,xf,"*linear");
+%! spl=interp1(xp,yp,xf,"*spline");
+%! cub=interp1(xp,yp,xf,"*cubic");
+%! near=interp1(xp,yp,xf,"*nearest");
+%! plot(xf,yf,";*original;",xf,near,";*nearest;",xf,lin,";*linear;",...
+%!      xf,cub,";*cubic;",xf,spl,";*spline;",xp,yp,"*;;");
+%! %--------------------------------------------------------
+%! % confirm that interpolated function matches the original
+
+%!shared xp, yp, xi, style
+%! xp=0:5;      yp = sin(2*pi*xp/5);
+%! xi = sort([-1, max(xp)*rand(1,6), max(xp)+1]);
+
+%!test style = "nearest";
+%!assert (interp1(xp, yp, [min(xp)-1, max(xp)+1]), [NaN, NaN]);
+%!assert (interp1(xp,yp,xp,style), yp, 100*eps);
+%!assert (interp1(xp,yp,xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp,style), yp, 100*eps);
+%!assert (isempty(interp1(xp',yp',[],style)));
+%!assert (isempty(interp1(xp,yp,[],style)));
+%!assert (interp1(xp,[yp',yp'],xi(:),style),...
+%!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,[yp',yp'],xi,style),
+%!	  interp1(xp,[yp',yp'],xi,["*",style]));
+
+%!test style = "linear";
+%!assert (interp1(xp, yp, [-1, max(xp)+1]), [NaN, NaN]);
+%!assert (interp1(xp,yp,xp,style), yp, 100*eps);
+%!assert (interp1(xp,yp,xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp,style), yp, 100*eps);
+%!assert (isempty(interp1(xp',yp',[],style)));
+%!assert (isempty(interp1(xp,yp,[],style)));
+%!assert (interp1(xp,[yp',yp'],xi(:),style),...
+%!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,[yp',yp'],xi,style),
+%!	  interp1(xp,[yp',yp'],xi,["*",style]),100*eps);
+
+%!test style = "cubic";
+%!assert (interp1(xp, yp, [-1, max(xp)+1]), [NaN, NaN]);
+%!assert (interp1(xp,yp,xp,style), yp, 100*eps);
+%!assert (interp1(xp,yp,xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp,style), yp, 100*eps);
+%!assert (isempty(interp1(xp',yp',[],style)));
+%!assert (isempty(interp1(xp,yp,[],style)));
+%!assert (interp1(xp,[yp',yp'],xi(:),style),...
+%!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,[yp',yp'],xi,style),
+%!	  interp1(xp,[yp',yp'],xi,["*",style]),1000*eps);
+
+%!test style = "spline";
+%!assert (interp1(xp, yp, [-1, max(xp) + 1]), [NaN, NaN]);
+%!assert (interp1(xp,yp,xp,style), yp, 100*eps);
+%!assert (interp1(xp,yp,xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp',style), yp', 100*eps);
+%!assert (interp1(xp',yp',xp,style), yp, 100*eps);
+%!assert (isempty(interp1(xp',yp',[],style)));
+%!assert (isempty(interp1(xp,yp,[],style)));
+%!assert (interp1(xp,[yp',yp'],xi(:),style),...
+%!	  [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)]);
+%!assert (interp1(xp,[yp',yp'],xi,style),
+%!	  interp1(xp,[yp',yp'],xi,["*",style]),10*eps);
+
+%!# test linear extrapolation
+%!assert (interp1([1:5],[3:2:11],[0,6],"linear","extrap"), [1, 13], eps);
+%!assert (interp1(xp, yp, [-1, max(xp)+1],"linear",5), [5, 5]);
+
+%!error interp1
+%!error interp1(1:2,1:2,1,"bogus")
+
+%!error interp1(1,1,1, "nearest");
+%!assert (interp1(1:2,1:2,1.4,"nearest"),1);
+%!error interp1(1,1,1, "linear");
+%!assert (interp1(1:2,1:2,1.4,"linear"),1.4);
+%!error interp1(1:3,1:3,1, "cubic");
+%!assert (interp1(1:4,1:4,1.4,"cubic"),1.4);
+%!error interp1(1:2,1:2,1, "spline");
+%!assert (interp1(1:3,1:3,1.4,"spline"),1.4);
+
+%!error interp1(1,1,1, "*nearest");
+%!assert (interp1(1:2:4,1:2:4,1.4,"*nearest"),1);
+%!error interp1(1,1,1, "*linear");
+%!assert (interp1(1:2:4,1:2:4,[0,1,1.4,3,4],"*linear"),[NaN,1,1.4,3,NaN]);
+%!error interp1(1:3,1:3,1, "*cubic");
+%!assert (interp1(1:2:8,1:2:8,1.4,"*cubic"),1.4);
+%!error interp1(1:2,1:2,1, "*spline");
+%!assert (interp1(1:2:6,1:2:6,1.4,"*spline"),1.4);
+
+%!assert (ppval(interp1(xp,yp,"nearest","pp"),xi),
+%!	  interp1(xp,yp,xi,"nearest","extrap"),10*eps);
+%!assert (ppval(interp1(xp,yp,"linear","pp"),xi),
+%!	  interp1(xp,yp,xi,"linear","extrap"),10*eps);
+%!assert (ppval(interp1(xp,yp,"cubic","pp"),xi),
+%!	  interp1(xp,yp,xi,"cubic","extrap"),10*eps);
+%!assert (ppval(interp1(xp,yp,"pchip","pp"),xi),
+%!	  interp1(xp,yp,xi,"pchip","extrap"),10*eps);
+%!assert (ppval(interp1(xp,yp,"spline","pp"),xi),
+%!	  interp1(xp,yp,xi,"spline","extrap"),10*eps);
+
+%!assert (ppval(interp1(xp,yp,"*nearest","pp"),xi),
+%!	  interp1(xp,yp,xi,"*nearest","extrap"),10*eps);
+%!assert (ppval(interp1(xp,yp,"*linear","pp"),xi),
+%!	  interp1(xp,yp,xi,"*linear","extrap"),10*eps);
+%!assert (ppval(interp1(xp,yp,"*cubic","pp"),xi),
+%!	  interp1(xp,yp,xi,"*cubic","extrap"),10*eps);
+%!assert (ppval(interp1(xp,yp,"*pchip","pp"),xi),
+%!	  interp1(xp,yp,xi,"*pchip","extrap"),10*eps);
+%!assert (ppval(interp1(xp,yp,"*spline","pp"),xi),
+%!	  interp1(xp,yp,xi,"*spline","extrap"),10*eps);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/interp2.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,316 @@
+## Copyright (C) 2000  Kai Habel
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{zi}=} interp2 (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
+## @deftypefnx {Function File} {@var{zi}=} interp2 (@var{Z}, @var{xi}, @var{yi})
+## @deftypefnx {Function File} {@var{zi}=} interp2 (@var{Z}, @var{n})
+## @deftypefnx {Function File} {@var{zi}=} interp2 (@dots{}, @var{method})
+## @deftypefnx {Function File} {@var{zi}=} interp2 (@dots{}, @var{method}, @var{extrapval})
+##
+## Two-dimensional interpolation. @var{x}, @var{y} and @var{z} describe a
+## surface function. If @var{x} and @var{y} are vectors their length
+## must correspondent to the size of @var{z}. @var{x} and @var{Yy must be
+## monotonic. If they are matrices they  must have the @code{meshgrid} 
+## format. 
+##
+## @table @code
+## @item interp2 (@var{x}, @var{y}, @var{Z}, @var{xi}, @var{yi}, @dots{}) 
+## Returns a matrix corresponding to the points described by the
+## matrices  @var{XI}, @var{YI}. 
+##
+## If the last argument is a string, the interpolation method can
+## be specified. The method can be 'linear', 'nearest' or 'cubic'.
+## If it is omitted 'linear' interpolation  is assumed.
+##
+## @item  interp2 (@var{z}, @var{xi}, @var{yi})
+## Assumes @code{@var{x} = 1:rows (@var{z})} and @code{@var{y} = 
+## 1:columns (@var{z})}
+## 
+## @item interp2 (@var{z}, @var{n}) 
+## Interleaves the Matrix @var{z} n-times. If @var{n} is ommited a value
+## of @code{@var{n} = 1} is assumed.
+## @end table
+##
+## The variable @var{method} defines the method to use for the
+## interpolation. It can take one of the values
+##
+## @table @asis
+## @item 'nearest'
+## Return the nearest neighbour.
+## @item 'linear'
+## Linear interpolation from nearest neighbours
+## @item 'pchip'
+## Piece-wise cubic hermite interpolating polynomial
+## @item 'cubic'
+## Cubic interpolation from four nearest neighbours
+## @item 'spline'
+## Cubic spline interpolation--smooth first and second derivatives
+## throughout the curve (Not implemented yet).
+## @end table
+##
+## If a scalar value @var{extrapval} is defined as the final value, then
+## values outside the mesh as set to this value. Note that in this case 
+## @var{method} must be defined as well. If @var{extrapval} is not
+## defined then NaN is assumed. 
+##
+## @seealso{interp1}
+## @end deftypefn
+
+## Author:	Kai Habel <kai.habel@gmx.de>
+## 2005-03-02 Thomas Weber <weber@num.uni-sb.de> 
+##     * Add test cases
+## 2005-03-02 Paul Kienzle <pkienzle@users.sf.net>
+##     * Simplify
+## 2005-04-23 Dmitri A. Sergatskov <dasergatskov@gmail.com>
+##     * Modified demo and test for new gnuplot interface
+## 2005-09-07 Hoxide <hoxide_dirac@yahoo.com.cn>
+##     * Add bicubic interpolation method
+##     * Fix the eat line bug when the last element of XI or YI is negative or zero.
+## 2005-11-26 Pierre Baldensperger <balden@libertysurf.fr>
+##     * Rather big modification (XI,YI no longer need to be
+##       "meshgridded") to be consistent with the help message
+##       above and for compatibility.
+
+
+function ZI = interp2 (varargin)
+  Z = X = Y = XI = YI = [];
+  n = [];
+  method = "linear";
+  extrapval = NaN;
+
+  switch nargin
+  case 1
+    Z = varargin{1};
+  case 2
+    if (ischar(varargin{2}))
+      [Z,method] = deal(varargin{:});
+    else
+      [Z,n] = deal(varargin{:});
+    endif
+  case 3
+    if (ischar(varargin{3}))
+      [Z,n,method] = deal(varargin{:});
+    else
+      [Z,XI,YI] = deal(varargin{:});
+    endif
+  case 4
+    if (ischar(varargin{4}))
+      [Z,XI,YI,method] = deal(varargin{:});
+    else
+      [Z,n,method,extrapval] = deal(varargin{:});
+    endif
+  case 5
+    if (ischar(varargin{4}))
+      [Z,XI,YI,method, extrapval] = deal(varargin{:});
+    else
+      [X,Y,Z,XI,YI] = deal(varargin{:});
+    endif
+  case 6 
+      [X,Y,Z,XI,YI,method] = deal(varargin{:});
+  case 7
+      [X,Y,Z,XI,YI,method,extrapval] = deal(varargin{:});
+  otherwise
+    print_usage ();
+  endswitch
+
+  ## Type checking.
+  if (!ismatrix(Z))
+    error("interp2 expected matrix Z"); 
+  endif
+  if (!isempty(n) && !isscalar(n))
+    error("interp2 expected scalar n"); 
+  endif
+  if (!ischar(method))
+    error("interp2 expected string 'method'"); 
+  endif
+  if (!isscalar(extrapval))
+    error("interp2 expected n extrapval");
+  endif
+
+  ## Define X,Y,XI,YI if needed
+  [zr, zc] = size (Z);
+  if (isempty(X))
+    X=[1:zc]; 
+    Y=[1:zr]; 
+  endif
+  if (!isnumeric(X) || !isnumeric(Y))
+    error("interp2 expected numeric X,Y"); 
+  endif
+  if (!isempty(n))
+    p=2^n; 
+    XI=[p:p*zc]/p; 
+    YI=[p:p*zr]'/p; 
+  endif
+  if (!isnumeric(XI) || !isnumeric(YI))
+    error("interp2 expected numeric XI,YI"); 
+  endif
+
+  ## If X and Y vectors produce a grid from them
+  if (isvector (X) && isvector (Y))
+    [X, Y] = meshgrid (X, Y);
+  elseif (! all(size (X) == size (Y)))
+    error ("X and Y must be matrices of same size");
+  endif
+  if (any(size (X) != size (Z)))
+    error ("X and Y size must match Z dimensions");
+  endif
+
+  ## If Xi and Yi are vectors of different orientation build a grid
+  if ((rows(XI)==1 && columns(YI)==1) || (columns(XI)==1 && rows(YI)==1))
+    [XI, YI] = meshgrid (XI, YI);
+  elseif (any(size(XI) != size(YI)))
+    error ("XI and YI must be matrices of same size");
+  endif
+
+  shape = size(XI);
+  XI = reshape(XI, 1, prod(shape));
+  YI = reshape(YI, 1, prod(shape));
+
+  xidx = lookup(X(1, 2:end-1), XI) + 1;
+  yidx = lookup(Y(2:end-1, 1), YI) + 1;
+
+  if (strcmp (method, "linear"))
+    ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy
+    ##
+    ## a-b
+    ## | |
+    ## c-d
+    a = Z(1:(zr - 1), 1:(zc - 1));
+    b = Z(1:(zr - 1), 2:zc) - a;
+    c = Z(2:zr, 1:(zc - 1)) - a;
+    d = Z(2:zr, 2:zc) - a - b - c;
+
+    idx = sub2ind(size(a),yidx,xidx);
+
+    ## scale XI,YI values to a 1-spaced grid
+    Xsc = (XI - X(1, xidx)) ./ (X(1, xidx + 1) - X(1, xidx));
+    Ysc = (YI - Y(yidx, 1)') ./ (Y(yidx + 1, 1) - Y(yidx, 1))';
+
+    ## apply plane equation
+    ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc;
+
+  elseif (strcmp (method, "nearest"))
+    xtable = X(1, :);
+    ytable = Y(:, 1)';
+    ii = (XI - xtable(xidx) > xtable(xidx + 1) - XI);
+    jj = (YI - ytable(yidx) > ytable(yidx + 1) - YI);
+    idx = sub2ind(size(Z),yidx+jj,xidx+ii);
+    ZI = Z(idx);
+
+  elseif (strcmp (method, "cubic"))
+    ## FIXME bicubic doesn't handle arbitrary XI, YI
+    ZI = bicubic(X, Y, Z, XI(1,:), YI(:,1));
+
+  elseif (strcmp (method, "spline"))
+    ## FIXME Implement 2-D (or in fact ND) spline interpolation
+    error ("interp2, spline interpolation not yet implemented");
+
+  else
+    error ("interpolation method not recognized");
+  endif
+
+  ## set points outside the table to NaN
+  ZI( XI < X(1,1) | XI > X(1,end) | YI < Y(1,1) | YI > Y(end,1) ) = extrapval;
+  ZI = reshape(ZI,shape);
+
+endfunction
+
+%!demo
+%! A=[13,-1,12;5,4,3;1,6,2];
+%! x=[0,1,4]; y=[10,11,12];
+%! xi=linspace(min(x),max(x),17);
+%! yi=linspace(min(y),max(y),26)';
+%! mesh(xi,yi,interp2(x,y,A,xi,yi,'linear'));
+%! [x,y] = meshgrid(x,y); 
+%! __gnuplot_raw__ ("set nohidden3d;\n")
+%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
+
+%!demo
+%! A=[13,-1,12;5,4,3;1,6,2];
+%! x=[0,1,4]; y=[10,11,12];
+%! xi=linspace(min(x),max(x),17);
+%! yi=linspace(min(y),max(y),26)';
+%! mesh(xi,yi,interp2(x,y,A,xi,yi,'nearest'));
+%! [x,y] = meshgrid(x,y); 
+%! __gnuplot_raw__ ("set nohidden3d;\n")
+%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
+
+%!#demo
+%! A=[13,-1,12;5,4,3;1,6,2];
+%! x=[0,1,2]; y=[10,11,12];
+%! xi=linspace(min(x),max(x),17);
+%! yi=linspace(min(y),max(y),26);
+%! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic'));
+%! [x,y] = meshgrid(x,y); 
+%! __gnuplot_raw__ ("set nohidden3d;\n")
+%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
+
+%!test % simple test
+%!  x = [1,2,3];
+%!  y = [4,5,6,7];
+%!  [X, Y] = meshgrid(x,y);
+%!  Orig = X.^2 + Y.^3;
+%!  xi = [1.2,2, 1.5];
+%!  yi = [6.2, 4.0, 5.0]';
+%!
+%!  Expected = ...
+%!    [243,   245.4,  243.9;
+%!      65.6,  68,     66.5;
+%!     126.6, 129,    127.5];
+%!  Result = interp2(x,y,Orig, xi, yi);
+%!
+%!  assert(Result, Expected, 1000*eps);
+
+%!test % 2^n form
+%!  x = [1,2,3];
+%!  y = [4,5,6,7];
+%!  [X, Y] = meshgrid(x,y);
+%!  Orig = X.^2 + Y.^3;
+%!  xi = [1:0.25:3]; yi = [4:0.25:7]';
+%!  Expected = interp2(x,y,Orig, xi, yi);
+%!  Result = interp2(Orig,2);
+%!  
+%!  assert(Result, Expected, 10*eps);
+
+%!test % matrix slice
+%!  A = eye(4);
+%!  assert(interp2(A,[1:4],[1:4]),[1,1,1,1]);
+
+%!test % non-gridded XI,YI
+%!  A = eye(4);
+%!  assert(interp2(A,[1,2;3,4],[1,3;2,4]),[1,0;0,1]);
+
+%!test % for values outside of boundaries
+%!  x = [1,2,3];
+%!  y = [4,5,6,7];
+%!  [X, Y] = meshgrid(x,y);
+%!  Orig = X.^2 + Y.^3;
+%!  xi = [0,4];
+%!  yi = [3,8]';
+%!  assert(interp2(x,y,Orig, xi, yi),[nan,nan;nan,nan]);
+%!  assert(interp2(x,y,Orig, xi, yi,'linear', 0),[0,0;0,0]);
+
+%!test % for values at boundaries
+%!  A=[1,2;3,4];
+%!  x=[0,1]; 
+%!  y=[2,3]';
+%!  assert(interp2(x,y,A,x,y,'linear'), A);
+%!  assert(interp2(x,y,A,x,y,'nearest'), A);
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/interpft.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,116 @@
+## Copyright (C) 2001 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} interpft (@var{x}, @var{n})
+## @deftypefnx {Function File} {} interpft (@var{x}, @var{n}, @var{dim})
+##
+## Fourier interpolation. If @var{x} is a vector, then @var{x} is
+## resampled with @var{n} points. The data in @var{x} is assumed to be
+## equispaced. If @var{x} is an array, then operate along each column of
+## the array seperately. If @var{dim} is specified, then interpolate
+## along the dimension @var{dim}.
+##
+## @code{interpft} assumes that the interpolated function is periodic,
+## and so assumption are made about the end points of the inetrpolation.
+##
+## @seealso{interp1}
+## @end deftypefn
+
+## Author: Paul Kienzle
+## 2001-02-11
+##    * initial version
+## 2002-03-17 aadler
+##    * added code to work on matrices as well 
+## 2006-05-25 dbateman
+##    * Make it matlab compatiable, cutting out the 2-D interpolation
+
+function z = interpft (x, n, dim)
+
+  if (nargin < 2 || nargin > 3)
+    print_usage ();
+  endif
+
+  if (nargin == 2)
+    if (isvector(x) && size(x,1) == 1)
+      dim = 2;
+    else
+      dim = 1;
+    endif
+  endif
+
+  if (!isscalar (n))
+    error ("interpft: n must be an integer scalar");
+  endif
+
+  nd = ndims(x);
+
+  if (dim < 1 || dim > nd)
+    error ("interpft: integrating over invalid dimension");
+  endif
+
+  perm = [dim:nd,1:(dim-1)];
+  x = permute(x, perm);
+  m = size (x, 1);
+
+  inc = 1;
+  while (inc * n < m)
+    inc++;
+  endwhile
+  y = fft (x) / m;
+  k = floor (m / 2);
+  sz = size(x);
+  sz(1) = n * inc - m;
+  idx = cell(nd,1);
+  for i = 2:nd
+    idx{i} = 1:sz(i);
+  endfor
+  idx{1} = 1:k;
+  z = cat (1, y(idx{:}), zeros(sz));
+  idx{1} = k+1:m;
+  z = cat (1, z, y(idx{:}));
+  z = n * ifft (z);
+
+  if (inc != 1)
+    sz(1) = n;
+    z = inc * reshape(z(1:inc:end),sz);
+  endif
+
+  z = ipermute (z, perm);
+endfunction
+
+%!demo
+%! t = 0 : 0.3 : pi; dt = t(2)-t(1);
+%! n = length (t); k = 100;
+%! ti = t(1) + [0 : k-1]*dt*n/k;
+%! y = sin (4*t + 0.3) .* cos (3*t - 0.1);
+%! yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1);
+%! plot (ti, yp, 'g;sin(4t+0.3)cos(3t-0.1);', ...
+%!       ti, interp1 (t, y, ti, 'spline'), 'b;spline;', ...
+%!       ti, interpft (y ,k), 'c;interpft;', ...
+%!       t, y, 'r+;data;');
+
+%!shared n,y
+%! x = [0:10]'; y = sin(x); n = length (x);
+%!assert (interpft(y, n), y, eps);
+%!assert (interpft(y', n), y', eps);
+%!assert (interpft([y,y],n), [y,y], eps);
+
+%!error (interpft(y,n,0))
+%!error (interpft(y,[n,n]))

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/polyarea.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,71 @@
+## Copyright (C) 1999 David M. Doolin <doolin@ce.berkeley.edu>
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} polyarea (@var{x}, @var{y})
+## @deftypefnx {Function File} {} polyarea (@var{x}, @var{y}, @var{dim})
+##
+## Determines area of a polygon by triangle method. The variables
+## @var{x} and @var{y} define the vertex pairs, and must therefore have
+## the same shape. Then might be either vectors or arrays. If they are
+## arrays then the columns of @var{x} and @var{y} are treated seperately
+## and an area returned for each.
+##
+## If the optional @var{dim} argument is given, then @code{polyarea}
+## works along this dimension of the arrays @var{x} and @var{y}.
+##
+## @end deftypefn
+
+## todo:  Add moments for centroid, etc.
+##
+## bugs and limitations:  
+##        Probably ought to be an optional check to make sure that
+##        traversing the vertices doesn't make any sides cross 
+##        (Is simple closed curve the technical definition of this?). 
+
+## Author: David M. Doolin <doolin@ce.berkeley.edu>
+## Date: 1999/04/14 15:52:20 $
+## Modified-by: 
+##    2000-01-15 Paul Kienzle <pkienzle@kienzle.powernet.co.uk>
+##    * use matlab compatible interface
+##    * return absolute value of area so traversal order doesn't matter
+##    2005-10-13 Torsten Finke
+##    * optimization saving half the sums and multiplies
+
+function a = polyarea (x, y, dim)
+  if (nargin != 2 && nargin != 3)
+    print_usage ();
+  elseif any (size(x) != size(y))
+    error ("polyarea: x and y must have the same shape");
+  else
+    if (nargin == 2)
+      a = abs (sum (x .* (shift (y, -1) - shift (y, 1)))) / 2;
+    else
+      a = abs (sum (x .* (shift (y, -1, dim) - shift (y, 1, dim)), dim)) / 2;
+    endif
+  endif
+endfunction
+
+%!shared x, y
+%! x = [1;1;3;3;1];
+%! y = [1;3;3;1;1];
+%!assert (polyarea(x,y), 4, eps)
+%!assert (polyarea([x,x],[y,y]), [4,4], eps)
+%!assert (polyarea([x,x],[y,y],1), [4,4], eps)
+%!assert (polyarea([x,x]',[y,y]',2), [4;4], eps)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/general/quadl.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,173 @@
+## Copyright (C) 1998 Walter Gautschi
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b})
+## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol})
+## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace})
+## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
+##
+## Numerically evaluate integral using adaptive Lobatto rule.
+## @code{quadl (@var{f}, @var{a}, @var{b})} approximates the integral of
+## @code{@var{f}(@var{x})} to machine precision. @var{f} is either a
+## function handle, inline function or string containing the name of
+## the function to evaluate. The function @var{f} must return a vector
+## of output values if given a vector of input values.
+##
+## If defined, @var{tol} defines the relative tolerance to which to
+## which to integrate @code{@var{f}(@var{x})}. While if @var{trace} is
+## defined, displays the left end point of the current interval, the 
+## interval length, and the partial integral.
+##
+## Additional arguments @var{p1}, etc, are passed directly to @var{f}.
+## To use default values for @var{tol} and @var{trace}, one may pass
+## empty matrices.
+##
+## Reference: W. Gander and W. Gautschi, 'Adaptive Quadrature - 
+## Revisited', BIT Vol. 40, No. 1, March 2000, pp. 84--101.
+## @url{http://www.inf.ethz.ch/personal/gander/}
+##
+## @end deftypefn
+
+##   Author: Walter Gautschi
+##   Date: 08/03/98
+##   Reference: Gander, Computermathematik, Birkhaeuser, 1992.
+
+## 2003-08-05 Shai Ayal
+##   * permission from author to release as GPL
+## 2004-02-10 Paul Kienzle
+##   * renamed to quadl for compatibility
+##   * replace global variable terminate2 with local function need_warning
+##   * add paper ref to docs
+
+function Q = quadl(f,a,b,tol,trace,varargin)
+  need_warning(1);
+  if (nargin < 4)
+    tol=[]; 
+  endif
+  if (nargin < 5)
+    trace = []; 
+  endif
+  if (isempty(tol))
+    tol = eps; 
+  endif
+  if (isempty(trace))
+    trace=0; 
+  endif
+  if (tol < eps)
+    tol = eps;
+  endif
+
+  m = (a+b)/2; 
+  h = (b-a)/2;
+  alpha = sqrt(2/3); 
+  beta = 1/sqrt(5);
+  x1 = .942882415695480; 
+  x2 = .641853342345781;
+  x3 = .236383199662150;
+  x = [a,m-x1*h,m-alpha*h,m-x2*h,m-beta*h,m-x3*h,m,m+x3*h,...
+       m+beta*h,m+x2*h,m+alpha*h,m+x1*h,b];
+  y = feval(f,x,varargin{:});
+  fa = y(1); 
+  fb = y(13);
+  i2 = (h/6)*(y(1)+y(13)+5*(y(5)+y(9)));
+  i1 = (h/1470)*(77*(y(1)+y(13))+432*(y(3)+y(11))+ ...
+		 625*(y(5)+y(9))+672*y(7));
+  is = h*(.0158271919734802*(y(1)+y(13))+.0942738402188500 ...
+	  *(y(2)+y(12))+.155071987336585*(y(3)+y(11))+ ...
+	  .188821573960182*(y(4)+y(10))+.199773405226859 ...
+	  *(y(5)+y(9))+.224926465333340*(y(6)+y(8))...
+	  +.242611071901408*y(7));
+  s = sign(is); 
+  if (s == 0),
+    s=1;
+  endif
+  erri1 = abs(i1-is);
+  erri2 = abs(i2-is);
+  R = 1; 
+  if (erri2 != 0)
+    R = erri1/erri2; 
+  endif
+  if (R > 0 && R < 1)
+    tol=tol/R; 
+  endif
+  is = s*abs(is)*tol/eps;
+  if (is == 0)
+    is = b-a;
+  endif
+  Q = adaptlobstp(f,a,b,fa,fb,is,trace,varargin{:});
+endfunction
+
+## ADAPTLOBSTP  Recursive function used by QUADL.
+##
+##   Q = ADAPTLOBSTP('F',A,B,FA,FB,IS,TRACE) tries to
+##   approximate the integral of F(X) from A to B to
+##   an appropriate relative error. The argument 'F' is
+##   a string containing the name of f.  The remaining
+##   arguments are generated by ADAPTLOB or by recursion.
+##
+##   Walter Gautschi, 08/03/98
+
+function Q = adaptlobstp(f,a,b,fa,fb,is,trace,varargin)
+  h = (b-a)/2; 
+  m = (a+b)/2;
+  alpha = sqrt(2/3); 
+  beta = 1/sqrt(5);
+  mll = m-alpha*h; 
+  ml = m-beta*h; 
+  mr = m+beta*h; 
+  mrr = m+alpha*h;
+  x = [mll,ml,m,mr,mrr];
+  y = feval(f,x,varargin{:});
+  fmll = y(1); 
+  fml = y(2); 
+  fm = y(3); 
+  fmr = y(4); 
+  fmrr = y(5);
+  i2 = (h/6)*(fa+fb+5*(fml+fmr));
+  i1 = (h/1470)*(77*(fa+fb)+432*(fmll+fmrr)+625*(fml+fmr) ...
+		 +672*fm);
+  if ((is+(i1-i2) == is) || (mll <= a) || (b <= mrr))
+    if (((m <= a) || (b <= m)) && need_warning())
+      warning(['Interval contains no more machine number. ',...
+               'Required tolerance may not be met.']);
+      need_warning(0);
+    endif
+    Q = i1;
+    if (trace)
+      disp([a b-a Q]);
+    endif
+  else
+    Q = adaptlobstp(f,a,mll,fa,fmll,is,trace,varargin{:})+...
+	adaptlobstp(f,mll,ml,fmll,fml,is,trace,varargin{:})+...
+	adaptlobstp(f,ml,m,fml,fm,is,trace,varargin{:})+...
+	adaptlobstp(f,m,mr,fm,fmr,is,trace,varargin{:})+...
+	adaptlobstp(f,mr,mrr,fmr,fmrr,is,trace,varargin{:})+...
+	adaptlobstp(f,mrr,b,fmrr,fb,is,trace,varargin{:});
+  endif
+endfunction
+
+function r = need_warning(v)
+  persistent w = [];
+  if (nargin == 0)
+    r = w;
+  else 
+    w = v; 
+  endif
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/miscellaneous/inputname.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,31 @@
+## Copyright (C) 2004 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## This program is free software and is in the public domain
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} inputname (@var{n})
+## 
+##    Return the text defining n-th input to the function.
+##
+## @end deftypefn
+
+function s = inputname (n)
+  s = evalin ('caller', sprintf('deblank(argn(%d,:));', n));
+endfunction
+
+## Warning: heap big magic in the following tests!!!
+## The test function builds a private context for each
+## test, with only the specified values shared between
+## them.  It does this using the following template:
+##
+##     function [<shared>] = testfn(<shared>)
+##        <test>
+##
+## To test inputname, I need a function context invoked
+## with known parameter names.  So define a couple of
+## shared parameters, et voila!, the test is trivial.
+%!shared hello,worldly
+%!assert(inputname(1),'hello');
+%!assert(inputname(2),'worldly');

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/__plt3__.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,65 @@
+## Copyright (C) 1996 John W. Eaton
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} __plt3__ (@var{x}, @var{y}, @var{z}, @var{fmt})
+## @end deftypefn
+
+## Author: Paul Kienzle <kienzle.powernet.co.uk>
+## 2001-04-06 Paul Kienzle <kienzle.powernet.co.uk>
+##     * __gnuplot_set__ nohidden3d; vector X,Y, matrix Z => meshgrid(X,Y)
+
+## Modified to use new gnuplot interface in octave > 2.9.0
+## Dmitri A. Sergatskov <dasergatskov@gmail.com>
+## April 18, 2005
+## Modified to use NaN as seperator for gnuplot, so multiple calls
+## aren't needed.
+## David Bateman <dbateman@free.fr>
+## May 25, 2006
+
+function __plt3__ (x, y, z, fmt)
+
+  if (isvector(x) && isvector(y))
+    if (isvector(z))
+      x = x(:); y=y(:); z=z(:);
+    elseif (length(x) == rows(z) && length(y) == columns(z))
+      error("plot3: [length(x), length(y)] must match size(z)");
+    else
+      [x,y] = meshgrid(x,y);
+    endif
+  endif
+
+  if (any(size(x) != size(y)) || any(size(x) != size(z)))
+    error("plot3: x, y, and z must have the same shape");
+  endif
+
+  unwind_protect
+    __gnuplot_set__  parametric;
+    __gnuplot_raw__ ("set nohidden3d;\n");
+
+    tmp = [([x; NaN*ones(1,size(x,2))])(:), ...
+	   ([y; NaN*ones(1,size(y,2))])(:), ...
+	   ([z; NaN*ones(1,size(z,2))])(:)];
+
+    cmd =  ["__gnuplot_splot__ tmp ", fmt, ";\n"];
+    eval (cmd);
+  unwind_protect_cleanup
+    __gnuplot_set__ noparametric; 
+  end_unwind_protect
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/ndgrid.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,74 @@
+## Copyright (C) 2006, Alexander Barth
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{y1}, @var{y2}, ...,  @var{y}n]} = ndgrid (@var{x1}, @var{x2}, ..., @var{x}n)
+## @deftypefnx {Function File} {[@var{y1}, @var{y2}, ...,  @var{y}n]} = ndgrid (@var{x})
+## Given n vectors @var{x1}, ... @var{x}n, ndgrid returns n arrays of dimension n. 
+## The elements of the ith output argument contains the elements of the vector 
+## @var{x}i repeated over all dimensions different from the ith dimension.
+## Calling ndgrid with only one input argument @var{x} is equivalent of calling ndgrid with 
+## all n input arguments equal to @var{x}:
+##
+## [@var{y1}, @var{y2}, ...,  @var{y}n] = ndgrid (@var{x}, ..., @var{x})
+## @seealso{meshgrid}
+## @end deftypefn
+
+## Author: Alexander Barth <abarth@marine.usf.edu>
+
+function varargout = ndgrid(varargin)
+
+  if (nargin == 1)
+    n = max([nargout 2]);  
+    ## If only one input argument is given, repeat it n-times
+    varargin{1:n} = varargin{1};
+  elseif (nargin >= nargout)
+    n = max([nargin 2]);  
+  else
+    error ("ndgrid: wrong number of input arguments");
+  endif
+
+
+  ## Determine the size of the output arguments
+  
+  shape = zeros(1,n);
+
+  for i=1:n
+    if (~isvector(varargin{i}))
+      error ("ndgrid: arguments must be vectors");
+    endif
+
+    shape(i) = length(varargin{i});
+  endfor
+
+  for i=1:n
+    ## size for reshape
+    r = ones(1,n);
+    r(i) = shape(i);
+
+    ## size for repmat
+    s = shape;
+    s(i) = 1;
+
+    varargout{i} = repmat(reshape(varargin{i},r) ,s);
+  endfor
+
+
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/plot/plot3.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,199 @@
+## Copyright (C) 1996 John W. Eaton
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} plot (@var{args})
+##
+## This function produces three-dimensional plots.  Many different
+## combinations of arguments are possible.  The simplest form is
+##
+## @example
+## plot3 (@var{x}, @var{y}, @var{z})
+## @end example
+##
+## @noindent
+## where the arguments are taken to be the vertices of the points to be
+## plotted in three dimensions. If all arguments are vectors of the same
+## length, then a single continuous line is drawn. If all arguments are
+## matrices, then each column of the matrices is treated as a seperate
+## line. No attempt is made to transpose the arguments to make the
+## number of rows match.
+##
+## To save a plot, in one of several image formats such as PostScript
+## or PNG, use the @code{print} command.
+##
+## An optional format argument can be given as
+##
+## @example
+## plot3 (@var{x}, @var{y}, @var{y}, @var{fmt})
+## @end example
+##
+## If the @var{fmt} argument is supplied, it is interpreted as
+## follows.  If @var{fmt} is missing, the default gnuplot line style
+## is assumed.
+##
+## @table @samp
+## @item -
+## Set lines plot style (default).
+##
+## @item .
+## Set dots plot style.
+##
+## @item @@
+## Set points plot style.
+##
+## @item -@@
+## Set linespoints plot style.
+##
+## @item ^
+## Set impulses plot style.
+##
+## @item L
+## Set steps plot style.
+##
+## @item @var{n}
+## Interpreted as the plot color if @var{n} is an integer in the range 1 to
+## 6.
+##
+## @item @var{nm}
+## If @var{nm} is a two digit integer and @var{m} is an integer in the
+## range 1 to 6, @var{m} is interpreted as the point style.  This is only
+## valid in combination with the @code{@@} or @code{-@@} specifiers.
+##
+## @item @var{c}
+## If @var{c} is one of @code{"k"}, @code{"r"}, @code{"g"}, @code{"b"},
+## @code{"m"}, @code{"c"}, or @code{"w"}, it is interpreted as the plot
+## color (black, red, green, blue, magenta, cyan, or white).
+##
+## @item ";title;"
+## Here @code{"title"} is the label for the key.
+##
+## @item +
+## @itemx *
+## @itemx o
+## @itemx x
+## Used in combination with the points or linespoints styles, set the point
+## style.
+## @end table
+##
+## The color line styles have the following meanings on terminals that
+## support color.
+##
+## @example
+## Number  Gnuplot colors  (lines)points style
+##   1       red                   *
+##   2       green                 +
+##   3       blue                  o
+##   4       magenta               x
+##   5       cyan                house
+##   6       brown            there exists
+## @end example
+##
+## The @var{fmt} argument can also be used to assign key titles.
+## To do so, include the desired title between semi-colons after the
+## formatting sequence described above, e.g. "+3;Key Title;"
+## Note that the last semi-colon is required and will generate an error if
+## it is left out.
+##
+## Arguments can also be given in groups of three as
+##
+## @example
+## plot3 (@var{x1}, @var{y1}, @var{y1}, @var{x2}, @var{y2}, @var{y2}, @dots{})
+## @end example
+## 
+## @noindent
+## where each set of three arguments are treated as seperate lines or
+## sets of lines in three dimensions.
+##
+## An example of the use of plot3 is
+##
+## @example
+## @group
+##    z = [0:0.05:5];
+##    plot3(cos(2*pi*z), sin(2*pi*z), z, ";helix;");
+## @end group
+## @end example
+##
+## @seealso{plot, semilogx, semilogy, loglog, polar, mesh, contour, __pltopt__
+## bar, stairs, errorbar, replot, xlabel, ylabel, title, print}
+## @end deftypefn
+
+## Author: Paul Kienzle
+##         (modified from __plt__.m)
+
+function plot3(varargin)
+
+  hold_state = ishold ();
+  
+  unwind_protect
+
+    x_set = 0;
+    y_set = 0;
+    z_set = 0;
+    
+    ## Gather arguments, decode format, and plot lines.
+    for arg = 1:length(varargin)
+      new = varargin{arg};
+      
+      if (ischar (new))
+	if (! z_set)
+	  error ("plot3: needs x, y, z");
+	endif
+	fmt = __pltopt__ ("plot3", new);
+	__plt3__(x, y, z, fmt);
+	hold on;
+	x_set = 0;
+	y_set = 0;
+	z_set = 0;
+      elseif (!x_set)
+	x = new;
+	x_set = 1;
+      elseif (!y_set)
+	y = new;
+	y_set = 1;
+      elseif (!z_set)
+	z = new;
+	z_set = 1;
+      else
+	__plt3__ (x, y, z, "");
+	hold on;
+	x = new;
+	y_set = 0;
+	z_set = 0;
+      endif
+       
+    endfor
+    
+    ## Handle last plot.
+    
+    if  (z_set)
+      __plt3__ (x, y, z, "");
+    elseif (x_set)
+      error ("plot3: needs x, y, z");
+    endif
+    
+  unwind_protect_cleanup
+    
+    if (! hold_state)
+      hold off;
+    endif
+    
+  end_unwind_protect
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/polynomial/pchip.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,154 @@
+## Copyright (C) 2001,2002  Kai Habel
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{pp} = } pchip (@var{x}, @var{y})
+## @deftypefnx {Function File} {@var{yi} = } pchip (@var{x}, @var{y}, @var{xi})
+##
+## Piecewise Cubic Hermite interpolating polynomial. Called with two
+## arguments, the piece-wise polynomial @var{pp} is returned, that may
+## later be used with @code{ppval} to evaluate the polynomial at
+## specific points.
+##
+## The variable @var{x} must be a strictly monotonic vector (either
+## increasing or decreasing). While @var{y} can be either a vector or
+## array. In the case where @var{y} is a vector, it must have a length
+## of @var{n}. If @var{y} is an array, then the size of @var{y} must
+## have the form
+## @iftex
+## @tex
+## $$[s_1, s_2, \cdots, s_k, n]$$
+## @end tex
+## @end iftex
+## @ifinfo
+## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]}
+## @end ifinfo
+## The array is then reshaped internally to a matrix where to leading
+## dimension is given by 
+## @iftex
+## @tex
+## $$s_1 s_2 \cdots s_k$$
+## @end tex
+## @end iftex
+## @ifinfo
+## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}}
+## @end ifinfo
+## and each row this matrix is then treated seperately. Note that this
+## is exactly the opposite treatment than @code{interp1} and is done
+## for compatiability.
+##
+## Called with a third input argument, @code{pchip} evaluates the 
+## piece-wise polynomial at the points @var{xi}. There is an equivalence
+## between @code{ppval (pchip (@var{x}, @var{y}), @var{xi})} and
+## @code{pchip (@var{x}, @var{y}, @var{xi})}.
+##
+## @seealso{spline, ppval, mkpp, unmkpp}
+## @end deftypefn
+
+## Author:  Kai Habel <kai.habel@gmx.de>
+## Date: 9. mar 2001
+##
+## S_k = a_k + b_k*x + c_k*x^2 + d_k*x^3; (spline polynom)
+##
+## 4 conditions:
+## S_k(x_k) = y_k;
+## S_k(x_k+1) = y_k+1;
+## S_k'(x_k) = y_k';
+## S_k'(x_k+1) = y_k+1';
+
+function ret = pchip (x, y, xi)
+
+  if (nargin < 2 || nargin > 3)
+    print_usage ();
+  endif
+
+  x = x(:);
+  n = length (x);
+
+  ## Check the size and shape of y
+  ndy = ndims (y);
+  szy = size (y);
+  if (ndy == 2 && (szy(1) == 1 || szy(2) == 1))
+    if (szy(1) == 1)
+      y = y';
+    else
+      szy = fliplr (szy);
+    endif
+  else
+    y = reshape (y, [prod(szy(1:end-1)), szy(end)])';
+  endif
+
+  h = diff(x);
+  if all(h<0)
+    x = flipud(x);
+    h = diff(x);
+    y = flipud(y);
+  elseif (any(h <= 0))
+    error("pchip: x must be strictly monotonic")
+  endif
+
+  if (rows(y) != n)
+    error("pchip: size of x and y must match");
+  endif
+
+  [ry, cy] = size (y);
+  if (cy > 1)
+    h = kron (diff (x), ones (1, cy));
+  endif
+  
+  dy = diff (y) ./ h;
+
+  a = y;
+  b = __pchip_deriv__(x,y);
+  c = - (b(2:n, :) + 2 * b(1:n - 1, :)) ./ h + 3 * diff (a) ./ h .^ 2;
+  d = (b(1:n - 1, :) + b(2:n, :)) ./ h.^2 - 2 * diff (a) ./ h.^3;
+
+  d = d(1:n - 1, :); c = c(1:n - 1, :);
+  b = b(1:n - 1, :); a = a(1:n - 1, :);
+  coeffs = [d(:), c(:), b(:), a(:)];
+  pp = mkpp (x, coeffs, szy(1:end-1));
+
+  if (nargin == 2)
+    ret = pp;
+  else
+    ret = ppval (pp, xi);
+  endif
+
+endfunction
+
+%!demo
+%! x = 0:8; 
+%! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0];
+%! xi = 0:0.01:8; 
+%! yspline = spline(x,y,xi);
+%! ypchip = pchip(x,y,xi);
+%! title("pchip and spline fit to discontinuous function");
+%! plot(xi,yspline,";spline;",...
+%!      xi,ypchip,"-;pchip;",...
+%!      x,y,"+;data;");
+%! %-------------------------------------------------------------------
+%! % confirm that pchip agreed better to discontinuous data than spline
+
+%!shared x,y
+%! x = 0:8; 
+%! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0];
+%!assert (pchip(x,y,x), y);
+%!assert (pchip(x,y,x'), y');
+%!assert (pchip(x',y',x'), y');
+%!assert (pchip(x',y',x), y);
+%!assert (isempty(pchip(x',y',[])));
+%!assert (isempty(pchip(x,y,[])));
+%!assert (pchip(x,[y;y],x), [pchip(x,y,x);pchip(x,y,x)])

--- a/scripts/polynomial/spline.m	Thu Jun 01 16:16:01 2006 +0000
+++ b/scripts/polynomial/spline.m	Thu Jun 01 19:05:32 2006 +0000
@@ -23,12 +23,12 @@
 ## @deftypefnx {Function File} {@var{yi} = } spline (@var{x}, @var{y}, @var{xi})
 ##
 ## Returns the cubic spline interpolation of @var{y} at the point
-## @var{x}. Called with two arguments the piece-wse polynomial @var{pp}
+## @var{x}. Called with two arguments the piece-wise polynomial @var{pp}
 ## that may later be used with @code{ppval} to evaluate the polynomial
 ## at specific points.
 ##
 ## The variable @var{x} must be a vector of length @var{n}, and @var{y}
-## can be either a vector of array. In the case where @var{y} is a
+## can be either a vector or array. In the case where @var{y} is a
 ## vector, it can have a length of either @var{n} or @code{@var{n} + 2}.
 ## If the length of @var{y} is @var{n}, then the 'not-a-knot' end
 ## condition is used. If the length of @var{y} is @code{@var{n} + 2},

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/sparse/pcg.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,482 @@
+## Copyright (C) 2004 Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{x} =} pcg (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M}, @var{x0}, @dots{})
+## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}, @var{eigest}] =} pcg (@dots{})
+##
+## Solves the linear system of equations @code{@var{A} * @var{x} =
+## @var{b}} by means of the  Preconditioned Conjugate Gradient iterative
+## method. The input arguments are
+##
+## @itemize
+## @item
+## @var{A} can be either a square (preferably sparse) matrix or a
+## function handle, inline function or string containing the name
+## of a function which computes @code{@var{A} * @var{x}}. In principle
+## @var{A} should be symmetric and positive definite; if @code{pcg}
+## finds @var{A} to not be positive definite, you will get a warning
+## message and the @var{flag} output parameter will be set.
+## 
+## @item
+## @var{b} is the right hand side vector.
+## 
+## @item
+## @var{tol} is the required relative tolerance for the residual error,
+## @code{@var{b} - @var{A} * @var{x}}. The iteration stops if @code{norm
+## (@var{b} - @var{A} * @var{x}) <= @var{tol} * norm (@var{b} - @var{A} *
+## @var{x0})}. If @var{tol} is empty or is omitted, the function sets
+## @code{@var{tol} = 1e-6} by default.
+## 
+## @item
+## @var{maxit} is the maximum allowable number of iterations; if
+## @code{[]} is supplied for @code{maxit}, or @code{pcg} has less
+## arguments, a default value equal to 20 is used.
+## 
+## @item
+## @var{M} is the (left) preconditioning matrix, so that the iteration is
+## (theoretically) equivalent to solving by @code{pcg} @code{@var{P} *
+## @var{x} = @var{M} \ @var{b}}, with @code{@var{P} = @var{M} \ @var{A}}.
+## Note that a proper choice of the preconditioner may dramatically
+## improve the overall performance of the method. Instead of matrix
+## @var{M}, the user may pass a function which returns the results of 
+## applying the inverse of @var{M} to a vector (usually this is the
+## preferred way of using the preconditioner). If @code{[]} is supplied
+## for @var{M}, or @var{M} is omitted, no preconditioning is applied.
+## 
+## @item
+## @var{x0} is the initial guess. If @var{x0} is empty or omitted, the 
+## function sets @var{x0} to a zero vector by default.
+## @end itemize
+## 
+## The arguments which follow @var{x0} are treated as parameters, and
+## passed in a proper way to any of the functions (@var{A} or @var{M})
+## which are passed to @code{pcg}. See the examples below for further
+## details. The output arguments are
+##
+## @itemize
+## @item
+## @var{x} is the computed approximation to the solution of
+## @code{@var{A} * @var{x} = @var{b}}.
+## 
+## @item
+## @var{flag} reports on the convergence. @code{@var{flag} = 0} means
+## the solution converged and the tolerance criterion given by @var{tol}
+## is satisfied. @code{@var{flag} = 1} means that the @var{maxit} limit
+## for the iteration count was reached. @code{@var{flag} = 3} reports that
+## the (preconditioned) matrix was found not positive definite.
+## 
+## @item
+## @var{relres} is the ratio of the final residual to its initial value,
+## measured in the Euclidean norm.
+## 
+## @item
+## @var{iter} is the actual number of iterations performed.
+##
+## @item 
+## @var{resvec} describes the convergence history of the method.
+## @code{@var{resvec} (i,1)} is the Euclidean norm of the residual, and
+## @code{@var{resvec} (i,2)} is the preconditioned residual norm,
+## after the (@var{i}-1)-th iteration, @code{@var{i} =
+## 1,2,...@var{iter}+1}. The preconditioned residual norm is defined as
+## @code{norm (@var{r}) ^ 2 = @var{r}' * (@var{M} \ @var{r})} where
+## @code{@var{r} = @var{b} - @var{A} * @var{x}}, see also the
+## description of @var{M}. If @var{eigest} is not required, only
+## @code{@var{resvec} (:,1)} is returned.
+## 
+## @item
+## @var{eigest} returns the estimate for the smallest @code{@var{eigest}
+## (1)} and largest @code{@var{eigest} (2)} eigenvalues of the
+## preconditioned matrix @code{@var{P} = @var{M} \ @var{A}}. In 
+## particular, if no preconditioning is used, the extimates for the
+## extreme eigenvalues of @var{A} are returned. @code{@var{eigest} (1)}
+## is an overestimate and @code{@var{eigest} (2)} is an underestimate, 
+## so that @code{@var{eigest} (2) / @var{eigest} (1)} is a lower bound
+## for @code{cond (@var{P}, 2)}, which nevertheless in the limit should
+## theoretically be equal to the actual value of the condition number. 
+## The method which computes @var{eigest} works only for symmetric positive
+## definite @var{A} and @var{M}, and the user is responsible for
+## verifying this assumption. 
+## @end itemize
+## 
+## Let us consider a trivial problem with a diagonal matrix (we exploit the
+## sparsity of A) 
+## 
+## @example
+## @group
+## 	N = 10; 
+## 	A = diag([1:N]); A = sparse(A);  
+## 	b = rand(N,1);
+## @end group
+## @end example
+## 
+## @sc{Example 1:} Simplest use of @code{pcg}
+## 
+## @example
+##   x = pcg(A,b)
+## @end example
+## 
+## @sc{Example 2:} @code{pcg} with a function which computes
+## @code{@var{A} * @var{x}}
+## 
+## @example
+## @group
+##   function y = applyA(x) 
+##     y = [1:N]'.*x; 
+##   endfunction
+##
+##   x = pcg('applyA',b)
+## @end group
+## @end example
+## 
+## @sc{Example 3:} Preconditioned iteration, with full diagnostics. The
+## preconditioner (quite strange, because even the original matrix
+## @var{A} is trivial) is defined as a function
+## 
+## @example
+## @group
+##   function y = applyM(x)		
+##     K = floor(length(x)-2); 
+##     y = x; 
+##     y(1:K) = x(1:K)./[1:K]';	
+##   endfunction
+## 
+##   [x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],[],'applyM')
+##   semilogy([1:iter+1], resvec);
+## @end group
+## @end example
+## 
+## @sc{Example 4:} Finally, a preconditioner which depends on a
+## parameter @var{k}.
+## 
+## @example
+## @group
+##   function y = applyM(x, varargin)
+##   K = varargin@{1@}; 
+##   y = x; y(1:K) = x(1:K)./[1:K]';	 
+##   endfuntion
+## 
+##   [x, flag, relres, iter, resvec, eigest] = ...
+##        pcg(A,b,[],[],'applyM',[],3)
+## @end group
+## @end example
+## 
+## @sc{References}
+## 
+## 	[1] C.T.Kelley, 'Iterative methods for linear and nonlinear equations',
+## 	SIAM, 1995 (the base PCG algorithm) 
+## 	
+## 	[2] Y.Saad, 'Iterative methods for sparse linear systems', PWS 1996
+## 	(condition number estimate from PCG) Revised version of this book is
+## 	available online at http://www-users.cs.umn.edu/~saad/books.html
+## 
+##
+## @seealso{sparse, pcr}
+## @end deftypefn
+
+## REVISION HISTORY
+##
+## 2004-05-21, Piotr Krzyzanowski:
+##	Added 4 demos and 4 tests
+##
+## 2004-05-18, Piotr Krzyzanowski:
+##	Warnings use warning() function now
+##
+## 2004-04-29, Piotr Krzyzanowski:
+##	Added more warning messages when FLAG is not required
+##
+## 2004-04-28, Piotr Krzyzanowski:
+## 	When eigest is required, resvec returns both the Euclidean and the
+##	preconditioned residual norm convergence history
+##
+## 2004-04-20, Piotr Krzyzanowski: 
+## 	Corrected eigenvalue estimator. Changed the tridiagonal matrix
+##	eigenvalue solver to regular eig
+## 
+
+function [x, flag, relres, iter, resvec, eigest] = ...
+		pcg( A, b, tol, maxit, M, x0, varargin )
+
+  if ((nargin < 6) || isempty(x0))
+    x = zeros(size(b));
+  else
+    x = x0;
+  endif
+
+  if (nargin < 5)
+    M = [];
+  endif
+
+  if ((nargin < 4) || isempty(maxit))
+    maxit = min(size(b,1),20);
+  endif
+
+  maxit = maxit + 2;
+
+  if ((nargin < 3) || isempty(tol))
+    tol = 1e-6;
+  endif
+
+  preconditioned_residual_out = false;
+  if (nargout > 5)
+    T = zeros(maxit,maxit);
+    preconditioned_residual_out = true;
+  endif
+
+  matrix_positive_definite = true;	# assume A is positive definite
+
+  p = zeros(size(b)); 
+  oldtau = 1; 
+  if (isnumeric(A))			# is A a matrix?
+    r = b - A*x; 
+  else					# then A should be a function!
+    r = b - feval(A,x,varargin{:});
+  endif
+
+  resvec(1,1) = norm(r); 
+  alpha = 1;
+  iter = 2;
+
+  while ((resvec(iter-1,1) > tol*resvec(1,1)) && (iter < maxit))
+    if (isnumeric(M))		# is M a matrix?
+      if isempty(M)		# if M is empty, use no precond
+	z = r;
+      else			# otherwise, apply the precond
+	z = M \ r;
+      endif
+    else			# then M should be a function!
+      z = feval(M,r,varargin{:});
+    endif
+    tau = z'*r; 
+    resvec(iter-1,2) = sqrt(tau);
+    beta = tau/oldtau;
+    oldtau = tau;
+    p = z + beta*p;
+    if (isnumeric(A))		# is A a matrix?
+      w = A*p;
+    else			# then A should be a function!
+      w = feval(A,p,varargin{:});
+    endif
+    oldalpha = alpha; 		# needed only for eigest
+    alpha = tau/(p'*w);
+    if (alpha <= 0.0) # negative matrix?
+      matrix_positive_definite = false;
+    endif
+    x = x + alpha*p;
+    r = r - alpha*w;
+    if ((nargout > 5) && (iter > 2))
+      T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ...
+	  [1 sqrt(beta); sqrt(beta) beta]./oldalpha;
+      ## EVS = eig(T(2:iter-1,2:iter-1));
+      ## fprintf(stderr,"PCG condest: %g (iteration: %d)\n", max(EVS)/min(EVS),iter);
+    endif
+    resvec(iter,1) = norm(r);
+    iter = iter + 1;
+  endwhile
+
+  if (nargout > 5)
+    if (matrix_positive_definite )
+      if (iter > 3)
+	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));
+      else
+	eigest = [NaN NaN];
+	warning("PCG: eigenvalue estimate failed: iteration converged too fast.");
+      endif
+    else
+      eigest = [NaN NaN];
+    endif
+
+    ## apply the preconditioner once more and finish with the precond
+    ## residual
+    if (isnumeric(M))		# is M a matrix?
+      if isempty(M)		# if M is empty, use no precond
+	z = r;
+      else			# otherwise, apply the precond
+	z = M\r;
+      endif
+    else			# then M should be a function!
+      z = feval(M,r,varargin{:});
+    endif
+    resvec(iter-1,2) = sqrt(r'*z);
+  else
+    resvec = resvec(:,1);	
+  endif
+
+  flag = 0;
+  relres = resvec(iter-1,1)./resvec(1,1);
+  iter = iter - 2;
+  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);
+    endif
+  else
+    if (nargout < 2)
+      fprintf(stderr, "PCG: converged in %d iterations. ", iter);
+      fprintf(stderr, "The initial residual norm was reduced %g times.\n",...
+	      1.0/relres);
+    endif
+  endif
+
+  if (!matrix_positive_definite)
+    flag = 3;
+    if (nargout < 2)
+      warning("PCG: matrix not positive definite?\n");
+    endif
+  endif
+endfunction
+
+%!demo
+%!
+%!	# Simplest usage of pcg (see also 'help pcg')
+%!
+%!	N = 10; 
+%!	A = diag([1:N]); b = rand(N,1); y = A\b; #y is the true solution
+%!  	x = pcg(A,b);
+%!	printf('The solution relative error is %g\n', norm(x-y)/norm(y));
+%!
+%!	# You shouldn't be afraid if pcg issues some warning messages in this
+%!	# example: watch out in the second example, why it takes N iterations 
+%!	# of pcg to converge to (a very accurate, by the way) solution
+%!demo
+%!
+%!	# Full output from pcg, except for the eigenvalue estimates
+%!	# We use this output to plot the convergence history  
+%!
+%!	N = 10; 
+%!	A = diag([1:N]); b = rand(N,1); X = A\b; #X is the true solution
+%!  	[x, flag, relres, iter, resvec] = pcg(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;');
+%!demo
+%!
+%!	# Full output from pcg, including the eigenvalue estimates
+%!	# Hilbert matrix is extremely ill conditioned, so pcg WILL have problems
+%!
+%!	N = 10; 
+%!	A = hilb(N); b = rand(N,1); X = A\b; #X is the true solution
+%!  	[x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],200);
+%!	printf('The solution relative error is %g\n', norm(x-X)/norm(X));
+%!	printf('Condition number estimate is %g\n', eigest(2)/eigest(1));
+%!	printf('Actual condition number is   %g\n', cond(A));
+%!	title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
+%!	semilogy([0:iter],resvec,['o-g;absolute residual;';'+-r;absolute preconditioned residual;']);
+%!demo
+%!
+%!	# Full output from pcg, including the eigenvalue estimates
+%!	# We use the 1-D Laplacian matrix for A, and cond(A) = O(N^2)
+%!	# and that's the reasone we need some preconditioner; here we take
+%!	# a very simple and not powerful Jacobi preconditioner, 
+%!	# which is the diagonal of A
+%!
+%!	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 = rand(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, eigest] = pcg(A,b,[],maxit);
+%!	printf('System condition number estimate is %g\n',eigest(2)/eigest(1));
+%!	title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
+%!	semilogy([0:iter],resvec(:,1),'o-g;NO preconditioning: absolute residual;');
+%!
+%!	pause(1);
+%!	# Test Jacobi preconditioner: it will not help much!!!
+%!
+%!	M = diag(diag(A)); # Jacobi preconditioner
+%!  	[x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],maxit,M);
+%!	printf('JACOBI preconditioned system condition number estimate is %g\n',eigest(2)/eigest(1));
+%!	hold on;
+%!	semilogy([0:iter],resvec(:,1),'o-r;JACOBI preconditioner: absolute residual;');
+%!
+%!	pause(1);
+%!	# Test nonoverlapping block Jacobi preconditioner: it will help much!
+%!
+%!	M = zeros(N,N);k=4
+%!	for i=1:k:N # form 1-D Laplacian matrix
+%!		M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1);
+%!	endfor
+%!  	[x, flag, relres, iter, resvec, eigest] = pcg(A,b,[],maxit,M);
+%!	printf('BLOCK JACOBI preconditioned system condition number estimate is %g\n',eigest(2)/eigest(1));
+%!	semilogy([0:iter],resvec(:,1),'o-b;BLOCK JACOBI preconditioner: absolute residual;');
+%!	hold off;
+%!test
+%!
+%!	#solve small diagonal system
+%!
+%!	N = 10; 
+%!	A = diag([1:N]); b = rand(N,1); X = A\b; #X is the true solution
+%!  	[x, flag] = pcg(A,b,[],N+1);
+%!	assert(norm(x-X)/norm(X),0,1e-10);
+%!	assert(flag,0);
+%!
+%!test
+%!
+%!	#solve small indefinite diagonal system
+%!	#despite A is indefinite, the iteration continues and converges
+%!	#indefiniteness of A is detected
+%!
+%!	N = 10; 
+%!	A = diag([1:N].*(-ones(1,N).^2)); b = rand(N,1); X = A\b; #X is the true solution
+%!  	[x, flag] = pcg(A,b,[],N+1);
+%!	assert(norm(x-X)/norm(X),0,1e-10);
+%!	assert(flag,3);
+%!
+%!test
+%!
+%!	#solve tridiagonal system, do not converge in default 20 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, eigest] = pcg(A,b,1e-12);
+%!	assert(flag);
+%!	assert(relres>1.0);
+%!	assert(iter,20); #should perform max allowable default number of iterations
+%!
+%!test
+%!
+%!	#solve tridiagonal system with 'prefect' preconditioner
+%!	#converges in one iteration, so the eigest does not work
+%!	#and issues a warning
+%!
+%!	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, eigest] = pcg(A,b,[],[],A,b);
+%!	assert(norm(x-X)/norm(X),0,1e-6);
+%!	assert(flag,0);
+%!	assert(iter,1); #should converge in one iteration
+%!	assert(isnan(eigest),isnan([NaN NaN]));
+%!

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/sparse/pcr.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,437 @@
+## Copyright (C) 2004 Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{x} =} pcr (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M}, @var{x0}, @dots{})
+## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} pcr (@dots{})
+## 
+## Solves the linear system of equations @code{@var{A} * @var{x} =
+## @var{b}} by means of the  Preconditioned Conjugate Residuals iterative
+## method. The input arguments are
+##
+## @itemize
+## @item
+## @var{A} can be either a square (preferably sparse) matrix or a
+## function handle, inline function or string containing the name
+## of a function which computes @code{@var{A} * @var{x}}. In principle
+## @var{A} should be symmetric and non-singular; if @code{pcr}
+## finds @var{A} to be numerically singular, you will get a warning
+## message and the @var{flag} output parameter will be set.
+## 
+## @item
+## @var{b} is the right hand side vector.
+## 
+## @item
+## @var{tol} is the required relative tolerance for the residual error,
+## @code{@var{b} - @var{A} * @var{x}}. The iteration stops if @code{norm
+## (@var{b} - @var{A} * @var{x}) <= @var{tol} * norm (@var{b} - @var{A} *
+## @var{x0})}. If @var{tol} is empty or is omitted, the function sets
+## @code{@var{tol} = 1e-6} by default.
+## 
+## @item
+## @var{maxit} is the maximum allowable number of iterations; if
+## @code{[]} is supplied for @code{maxit}, or @code{pcr} has less
+## arguments, a default value equal to 20 is used.
+##
+## @item
+## @var{M} is the (left) preconditioning matrix, so that the iteration is
+## (theoretically) equivalent to solving by @code{pcr} @code{@var{P} *
+## @var{x} = @var{M} \ @var{b}}, with @code{@var{P} = @var{M} \ @var{A}}.
+## Note that a proper choice of the preconditioner may dramatically
+## improve the overall performance of the method. Instead of matrix
+## @var{M}, the user may pass a function which returns the results of 
+## applying the inverse of @var{M} to a vector (usually this is the
+## preferred way of using the preconditioner). If @code{[]} is supplied
+## for @var{M}, or @var{M} is omitted, no preconditioning is applied.
+## 
+## @item
+## @var{x0} is the initial guess. If @var{x0} is empty or omitted, the 
+## function sets @var{x0} to a zero vector by default.
+## @end itemize
+## 
+## The arguments which follow @var{x0} are treated as parameters, and
+## passed in a proper way to any of the functions (@var{A} or @var{M})
+## which are passed to @code{pcr}. See the examples below for further
+## details. The output arguments are
+##
+## @itemize
+## @item
+## @var{x} is the computed approximation to the solution of
+## @code{@var{A} * @var{x} = @var{b}}.
+## 
+## @item
+## @var{flag} reports on the convergence. @code{@var{flag} = 0} means
+## the solution converged and the tolerance criterion given by @var{tol}
+## is satisfied. @code{@var{flag} = 1} means that the @var{maxit} limit
+## for the iteration count was reached. @code{@var{flag} = 3} reports t
+## @code{pcr} breakdown, see [1] for details.
+## 
+## @item
+## @var{relres} is the ratio of the final residual to its initial value,
+## measured in the Euclidean norm.
+## 
+## @item
+## @var{iter} is the actual number of iterations performed.
+##
+## @item 
+## @var{resvec} describes the convergence history of the method,
+## so that @code{@var{resvec} (i)} contains the Euclidean norms of the 
+## residualafter the (@var{i}-1)-th iteration, @code{@var{i} =
+## 1,2,...@var{iter}+1}.
+## @end itemize
+## 
+## Let us consider a trivial problem with a diagonal matrix (we exploit the
+## sparsity of A) 
+## 
+## @example
+## @group
+## 	N = 10; 
+## 	A = diag([1:N]); A = sparse(A);  
+## 	b = rand(N,1);
+## @end group
+## @end example
+## 
+## @sc{Example 1:} Simplest use of @code{pcr}
+## 
+## @example
+##   x = pcr(A, b)
+## @end example
+## 
+## @sc{Example 2:} @code{pcr} with a function which computes
+## @code{@var{A} * @var{x}}.
+##
+## @example
+## @group
+##   function y = applyA(x) 
+##     y = [1:10]'.*x; 
+##   endfunction
+## 
+##   x = pcr('applyA',b)
+## @end group
+## @end example
+## 
+## @sc{Example 3:}  Preconditioned iteration, with full diagnostics. The
+## preconditioner (quite strange, because even the original matrix
+## @var{A} is trivial) is defined as a function
+## 
+## @example
+## @group
+##   function y = applyM(x)		
+##     K = floor(length(x)-2); 
+##     y = x; 
+##     y(1:K) = x(1:K)./[1:K]';	
+##   endfunction
+## 
+##   [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM')
+##   semilogy([1:iter+1], resvec);
+## @end group
+## @end example
+##
+## @sc{Example 4:} Finally, a preconditioner which depends on a
+## parameter @var{k}.
+## 
+## @example
+## @group
+##   function y = applyM(x, varargin)
+##     K = varargin@{1@}; 
+##     y = x; y(1:K) = x(1:K)./[1:K]';	 
+##   endfunction
+## 
+##   [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM',[],3)
+## @end group
+## @end example
+## 
+## @sc{References}
+## 
+## 	[1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of
+##  	Equations", section 9.5.4; Springer, 1994
+##
+## @seealso{sparse, pcg}
+## @end deftypefn
+
+## REVISION HISTORY
+##
+## 2004-08-14, Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
+## 
+## 	Added 4 demos and 4 tests
+##  
+## 2004-08-01, Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
+## 
+## 	First version of pcr code
+
+function [x, flag, relres, iter, resvec] = pcr(A,b,tol,maxit,M,x0,varargin)
+
+  breakdown = false;
+
+  if ((nargin < 6) || isempty(x0))
+    x = zeros(size(b));
+  else
+    x = x0;
+  endif
+
+  if (nargin < 5)
+    M = [];
+  endif
+
+  if ((nargin < 4) || isempty(maxit))
+    maxit = 20;
+  endif
+
+  maxit = maxit + 2;
+
+  if ((nargin < 3) || isempty(tol))
+    tol = 1e-6;
+  endif
+
+  if (nargin < 2)
+    print_usage();
+  endif
+
+  ##  init
+  if (isnumeric(A))		# is A a matrix?
+    r = b-A*x;
+  else				# then A should be a function!
+    r = b-feval(A,x,varargin{:});
+  endif
+
+  if (isnumeric(M))		# is M a matrix?
+    if isempty(M)		# if M is empty, use no precond
+      p = r;
+    else			# otherwise, apply the precond
+      p = M\r;
+    endif
+  else				# then M should be a function!
+    p = feval(M,r,varargin{:});
+  endif
+
+  iter = 2;
+
+  b_bot_old = 1;
+  q_old = p_old = s_old = zeros(size(x));
+
+  if (isnumeric(A))		# is A a matrix?
+    q = A*p;
+  else				# then A should be a function!
+    q = feval(A,p,varargin{:});
+  endif
+	
+  resvec(1) = abs(norm(r)); 
+
+  ## iteration
+  while ((resvec(iter-1) > tol*resvec(1)) && (iter < maxit))
+	
+    if (isnumeric(M))		# is M a matrix?
+      if isempty(M)		# if M is empty, use no precond
+	s = q;
+      else			# otherwise, apply the precond
+	s = M\q;
+      endif
+    else			# then M should be a function!
+      s = feval(M,q,varargin{:});
+    endif
+    b_top = r'*s;
+    b_bot = q'*s;
+	
+    if (b_bot == 0.0)
+      breakdown = true;
+      break;
+    endif
+    lambda = b_top/b_bot;
+	
+    x = x + lambda*p;
+    r = r - lambda*q;
+	
+    if (isnumeric(A))		# is A a matrix?
+      t = A*s;
+    else			# then A should be a function!
+      t = feval(A,s,varargin{:});
+    endif
+	
+    alpha0 = (t'*s)/b_bot;
+    alpha1 = (t'*s_old)/b_bot_old;
+	
+    p_temp = p; q_temp = q;
+    p = s - alpha0*p - alpha1*p_old;
+    q = t - alpha0*q - alpha1*q_old;
+	
+    s_old = s;
+    p_old = p_temp;
+    q_old = q_temp;
+    b_bot_old = b_bot;
+	
+	
+    resvec(iter) = abs(norm(r));
+    iter = iter + 1;
+  endwhile
+
+  flag = 0;
+  relres = resvec(iter-1)./resvec(1);
+  iter = iter - 2;
+  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);
+    endif
+  else
+    if ((nargout < 2) && (~breakdown))
+      fprintf(stderr, "PCR: converged in %d iterations. \n", iter);
+      fprintf(stderr, "The initial residual norm was reduced %g times.\n",...
+	      1.0/relres);
+    endif
+  endif
+  if (breakdown)
+    flag = 3;
+    if (nargout < 2)
+      warning("PCR: breakdown occured.\n");
+      warning("System matrix singular or preconditioner indefinite?\n");
+    endif
+  endif
+
+endfunction
+
+%!demo
+%!
+%!	# 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
+%!
+%!	# 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;');
+%!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
+%!
+%!	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;');
+%!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
+%!
+%!	# 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!
+%!
+%!  	[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!!!
+%!
+%!	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 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);
+%!
+%!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);
+%!
+%!test
+%!
+%!	#solve tridiagonal system with 'prefect' 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
+%!

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/strings/mat2str.m	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,113 @@
+## Copyright (C) 2002 Rolf Fabian <fabian@tu-cottbus.de>
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{s} =} mat2str (@var{x}, @var{n})
+## @deftypefnx {Function File} {@var{s} =} mat2str (@dots{}, 'class')
+##
+## Format real/complex numerial matrices as strings. This function
+## returns values that are suitable for the use of the @code{eval}
+## function.
+##
+## The precision of the values is given by @var{n}. If @var{n} is a
+## scalar then both real and imaginary parts of the matrix are printed
+## to the same precision. Otherwise @code{@var{n} (1)} defines the
+## precision of the real part and @code{@var{n} (2)} defines the
+## precision of the imaginary part. The default for @var{n} is 17.
+##
+## If the argument 'class' is given, then the class of @var{x} is
+## included in the string in such a way that the eval will result in the
+## construction of a matrix of the same class.
+##
+## @example
+## @group
+##    mat2str( [ -1/3 + i/7; 1/3 - i/7 ], [4 2] )
+##    @result{} '[-0.3333+0.14i;0.3333-0.14i]'
+##    mat2str( [ -1/3 +i/7; 1/3 -i/7 ], [4 2] )
+##    @result{} '[-0.3333+0i,0+0.14i;0.3333+0i,-0-0.14i]'
+##    mat2str( int16([1 -1]), 'class')
+##    @result{} 'int16([1,-1])'
+## @end group
+## @end example
+##
+## @seealso{sprintf, int2str}
+## @end deftypefn
+
+function s = mat2str(x,n,cls)
+
+  if (nargin < 2 || isempty(n))
+    n = 17;		   # default precision
+  endif
+
+  if (nargin < 3)
+    if (ischar(n))
+      cls = n;
+      n = 17;
+    else
+      cls = '';
+    endif
+  endif
+
+  if (nargin < 1 || nargin > 3 || ischar(x) || isstruct(x) ||
+     ischar(n) || isstruct(n) || isstruct(cls))
+    usage ("mat2str");
+  endif
+  
+  if (ndims (x) > 2)
+    error ("mat2str: x must be two dimensional");
+  endif
+
+  if (!(COMPLEX = is_complex(x)))
+    FMT = sprintf("%%.%dg", n(1));
+  else
+    if (length(n) == 1 )
+      n = [n, n];
+    endif
+    FMT = sprintf("%%.%dg%%+.%dgi", n(1), n(2));
+  endif
+
+  [nr, nc] = size(x);
+
+  if (nr*nc == 0)         # empty .. only print brackets
+    s = "[]";
+  elseif (nr*nc == 1)	# scalar x .. don't print brackets
+    if (!COMPLEX)
+      s = sprintf( FMT, x );
+    else
+      s = sprintf( FMT, real(x), imag(x) );
+    endif
+  else			# non-scalar x .. print brackets
+    FMT = [FMT, ','];
+    if (!COMPLEX)
+      s = sprintf( FMT, x.' );
+    else
+      x = x.';
+      s = sprintf( FMT, [ real(x(:))'; imag(x(:))' ] );
+    endif
+
+    s = ["[", s];
+    s (length(s)) = "]";
+    IND = find(s == ",");
+    s (IND(nc:nc:length(IND)) ) = ";";
+  endif
+
+  if (strcmp ("class", cls))
+    s = [class(x), "(", s, ")"]
+  endif
+endfunction

--- a/src/ChangeLog	Thu Jun 01 16:16:01 2006 +0000
+++ b/src/ChangeLog	Thu Jun 01 19:05:32 2006 +0000
@@ -1,3 +1,7 @@
+2006-06-01  David Bateman  <dbateman@free.fr>
+
+	* DLD-FUNCTIONS/__pchip_deriv__.cc: New file.
+
 2006-05-31  John W. Eaton  <jwe@octave.org>
 
 	* load-path.h (load_path::set_command_line_path): Make it

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/__pchip_deriv__.cc	Thu Jun 01 19:05:32 2006 +0000
@@ -0,0 +1,96 @@
+/* 
+
+Copyright (C) 2002  Kai Habel
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "oct-obj.h"
+#include "utils.h"
+#include "f77-fcn.h"
+
+extern "C" {
+  extern int F77_FUNC (dpchim, DPCHIM)
+    (const int &n, double *x, double *f, 
+	  double *d, const int &incfd, int *ierr);
+}
+
+DEFUN_DLD(__pchip_deriv__, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} __pchip_deriv__ (@var{x}, @var{y})\n\
+Wrapper for SLATEC/PCHIP function DPCHIM to calculate the derivates for\n\
+piecewise polynomials. You should be using @code{pchip} function instead.\n\
+@end deftypefn")
+{
+  octave_value retval;
+  const int nargin = args.length();
+
+  if (nargin == 2)
+    {
+      ColumnVector xvec(args(0).vector_value());
+      Matrix ymat(args(1).matrix_value());
+      int nx = xvec.length();
+      int nyr = ymat.rows();
+      int nyc = ymat.columns();
+
+      if (nx != nyr)
+        {
+          error("number of rows for x and y must match");
+          return retval;
+        }
+
+      ColumnVector dvec(nx),yvec(nx);
+      Matrix dmat(nyr,nyc);
+
+      int ierr;
+      const int incfd = 1;
+      for (int c = 0; c < nyc; c++)
+        {
+          for (int r=0; r<nx; r++)
+            {
+              yvec(r) = ymat(r,c);
+            }
+
+          F77_FUNC (dpchim, DPCHIM) (nx, xvec.fortran_vec(), 
+				     yvec.fortran_vec(), 
+				     dvec.fortran_vec(), incfd, &ierr);
+
+	  if (ierr < 0)
+            {
+	      error ("DPCHIM error: %i\n",ierr);
+              return retval;
+            }
+          for (int r=0; r<nx; r++)
+            {
+              dmat(r,c) = dvec(r);
+            }
+        }
+
+      retval = dmat;
+    }
+
+  return retval;
+}