function B = biharmonic_op(m,Lx,n,Ly)
%
% Returns the sparse matrix for the discrete biharmonic operator on a
% regular mxn grid of interior points on a rectangular domain
% [0,Lx]x[0,Ly], assuming Dirichlet boundary conditions and the
% standard treatment of the normal derivative boundary conditions,
% i.e., reflecting the first interior gridline across the boundary on
% each side.
%
% The inputs are as follows: 
%
% Required:
%   m            = number of meshpoints in the x-direction. 
%
% Optional: 
%   Lx          = length of rectangle in the x-direction (default 1). 
%   n           = number of meshpoints in the y-direction (default m).
%   Ly          = length of rectangle in the y-direction (default 1 or Lx). 
%
% The output is 
%   B  = sparse discrete biharmonic operator. 
%

% I am indebted to Professor Homer F. Walker,
% Mathematical Sciences Department, Worchester Polytechnic Institute,
% who provided this code in a private communication. 

if nargin > 4
   error('biharmonic_op: Number of arguments must be <= 4.');
end 

if nargin == 1,
  n = m; Lx = 1; Ly = 1;
elseif nargin == 2, 
  n = m; Ly = Lx;
elseif nargin == 3,
  Ly = Lx;
end
if isempty(Lx), Lx = 1; end
if isempty(n), n = m; end
if isempty(Ly), Ly = Lx; end

hx = (m+1)\Lx;  
hy = (n+1)\Ly;

Im = speye(m,m);
E1m = sparse(2:m,1:m-1,1,m,m);
E2m = sparse(3:m,1:m-2,1,m,m);

Mx = 6*Im -4*(E1m+E1m') + (E2m+E2m');
Mx(1,1) = Mx(1,1) + 1;
Mx(m,m) = Mx(m,m) + 1;

In = speye(n,n);
E1n = sparse(2:n,1:n-1,1,n,n);
E2n = sparse(3:n,1:n-2,1,n,n);

My = 6*In -4*(E1n+E1n') + (E2n+E2n');
My(1,1) = My(1,1) + 1;
My(n,n) = My(n,n) + 1;

Dm = 4*Im - 2*(E1m+E1m');
Dn = 2*In - (E1n+E1n');

B = ((hx^4)\1)*kron(In,Mx) + ((hy^4)\1)*kron(My,Im) + (((hx^2)*(hy^2))\1)*kron(Dn,Dm);

