--- a/doc/interpreter/interp.txi	Thu Jun 14 01:18:26 2007 +0000
+++ b/doc/interpreter/interp.txi	Thu Jun 14 06:56:42 2007 +0000
@@ -15,6 +15,43 @@
 
 @DOCSTRING(interp1)
 
+There are some important differences between the various interpolation
+methods. The 'spline' method enforces that both the first and second
+derivatives of the interpolated values have a continuous derivative,
+whereas the other methods do not. This can be demonstrated by the code
+
+@example
+@group
+t = 0 : 0.3 : pi; dt = t(2)-t(1);
+n = length (t); k = 100; dti = dt*n/k;
+ti = t(1) + [0 : k-1]*dti;
+y = sin (4*t + 0.3) .* cos (3*t - 0.1);
+ddyc = diff(diff(interp1(t,y,ti,'cubic'))./dti)./dti;
+ddys = diff(diff(interp1(t,y,ti,'spline'))./dti)./dti;
+ddyp = diff(diff(interp1(t,y,ti,'pchip'))./dti)./dti;
+plot (ti(2:end-1), ddyc,'g+',ti(2:end-1),ddys,'b*', ...
+      ti(2:end-1),ddyp,'c^');
+legend('cubic','spline','pchip');
+@end group
+@end example
+
+@ifnotinfo
+@noindent
+The result of which can be seen in @ref{fig:interpderiv}.
+
+@float Figure,fig:interpderiv
+@image{interpderiv,8cm}
+@caption{Comparison of second derivative of interpolated values for
+various interpolation methods}
+@end float
+@end ifnotinfo
+
+This means that in general the 'spline' method results in smooth
+results. If the function to be interpolated is in fact smooth, then
+'spline' will give excellent results. However, if the function to be
+evaluated is in some manner discontinuous, then 'cubic' or 'pchip'
+interpolation might give better results.
+
 Fourier interpolation, is a resampling technique where a signal is
 converted to the frequency domain, padded with zeros and then
 reconverted to the time domain.
@@ -40,8 +77,20 @@
 @end group
 @end example
 
+@ifinfo
 which demonstrates the poor behavior of Fourier interpolation for non
 periodic functions.
+@end ifinfo
+@ifnotinfo
+which demonstrates the poor behavior of Fourier interpolation for non
+periodic functions, as can be seen in @ref{fig:interpft}.
+
+@float Figure,fig:interpft
+@image{interpft,8cm}
+@caption{Comparison of @code{interp1} and @code{interpft} for non
+periodic data}
+@end float
+@end ifnotinfo
 
 In additional the support function @code{spline} and @code{lookup} that
 underlie the @code{interp1} function can be called directly.
@@ -81,11 +130,12 @@
 f = @@(x,y,z) x.^2 - y - z.^2;
 [xx, yy, zz] = meshgrid (x, y, z);
 v = f (xx,yy,zz);
-xi = yi = zi = -1:0.5:1;
+xi = yi = zi = -1:0.1:1;
 [xxi, yyi, zzi] = meshgrid (xi, yi, zi);
-vi = interp3(x, y, z, v, xxi, yyi, zzi);
+vi = interp3(x, y, z, v, xxi, yyi, zzi, 'spline');
 [xxi, yyi, zzi] = ndgrid (xi, yi, zi);
-vi2 = interpn(x, y, z, v, xxi, yyi, zzi);
+vi2 = interpn(x, y, z, v, xxi, yyi, zzi, 'spline');
+mesh (yi, zi, squeeze (vi2(1,:,:)));
 @end group
 @end example
 
@@ -93,6 +143,14 @@
 where @code{vi} and @code{vi2} are identical. The reversal of the
 dimensions is treated in the @code{meshgrid} and @code{ndgrid} functions
 respectively.
+@ifnotinfo
+The result of this code can be seen in @ref{fig:interpn}.
+
+@float Figure,fig:interpn
+@image{interpn,8cm}
+@caption{Demonstration of the use of @code{interpn}}
+@end float
+@end ifnotinfo
 
 In additional the support function @code{bicubic} that underlies the
 cubic interpolation of @code{interp2} function can be called directly.