From octave-sources-request at bevo dot che dot wisc dot edu  Tue Mar 30 20:34:14 1999
Subject: fd stencil utility
From: David Doolin <doolin at cs dot utk dot edu>
To: octave-sources at bevo dot che dot wisc dot edu
cc: doolin at cs dot utk dot edu
Date: Tue, 30 Mar 1999 21:34:10 -0500


Returns a tridiagonal fd matrix.  This definitely
falls into the category of "homework ware".  So far 
I have only done 1D problems, so there be a modification
to this coming pretty quick.

Feedback appreciated.


Dave Doolin

doolin at ce dot berkeley dot edu
doolin at cs dot utk dot edu

----------------cut here----------------

% [ retval ] = fdstencil(dim, numnodes)
%
% input:
%       dim: the dimension of the stencil 
%
% output:
%       retval:  the finite difference stencil
%
% bugs and limitations:  

% This script properly belongs in octave/scripts/special-matrix/
%
% This code is copyright David M. Doolin under the GNU
% General Public License (GPL). 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 file.  If not, write to the Free Software Foundation,
% 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
%
% $Author: doolin $  doolin at ce dot berkeley dot edu
% $Date: 1999/03/31 02:23:49 $
% $Source: /shag/homes/doolin/cvsroot/octave/finitediffs/fdstencil.m,v $
% $Revision: 1.2 $



function [ retval ] = fdstencil(dim, numnodes)

if (nargin != 2)
  usage("fdstencil(dim, numnodes)");
endif

if ( (dim <= 0) || (dim>3) )
  error("fdstencil: expecting argument 1 <= dim <= 3");
endif

retval = zeros(numnodes, numnodes);

for i = 1:dim

% Get the diagonal matrices.
    diag_m1 = ones(numnodes-1,1);
    diag_0 = -2*ones(numnodes,1);
    dm1 = diag(diag_m1, -1);
    d0 = diag(diag_0, 0);
    dp1 = diag(diag_m1,1);
    retval = retval + dm1+d0+dp1;
end %% for 

endfunction