function [L, B, mx, my] = plate_matrices(nx, ny, dx, dy, bc)
% PLATE_MATRICES  Sparse Laplacian and biharmonic matrices on the interior nodes.
%
%   [L, B, mx, my] = plate_matrices(nx, ny, dx, dy, bc) assembles, for the
%   mx = nx-2 by my = ny-2 interior nodes of the grid (edge nodes carry w = 0),
%     L   the five-point Laplacian            (5 nonzeros per row)
%     B   the thirteen-point biharmonic       (13 nonzeros per row)
%   as sparse matrices acting on w(:) with the interior field stored as an
%   my x mx array, so the y index runs fastest, matching w(2:end-1, 2:end-1)(:).
%
%   The stencils are the same as the ghost-node scheme in plate_fd: the second
%   difference with w = 0 at the edge is D2; the fourth difference is D2*D2 for
%   simply supported edges (odd reflection, w_nn = 0), and for clamped edges the
%   even reflection w_ghost = w_1 changes the end coefficient of the fourth
%   difference from 5/dx^4 to 7/dx^4.  The mixed term 2 w_xxyy never reaches a
%   ghost node.  A test in the lab checks B against the explicit stencil.
%
%   Everything is built from Kronecker products so the matrices stay sparse.

mx = nx - 2;
my = ny - 2;
ex = ones(mx, 1);
ey = ones(my, 1);
D2x = spdiags([ex -2 * ex ex], -1:1, mx, mx) / dx^2;
D2y = spdiags([ey -2 * ey ey], -1:1, my, my) / dy^2;
Ix = speye(mx);
Iy = speye(my);

D4x = D2x * D2x;
D4y = D2y * D2y;
switch bc
    case 'simply-supported'
        % D2*D2 already encodes the odd reflection
    case 'clamped'
        D4x(1, 1) = 7 / dx^4;   D4x(mx, mx) = 7 / dx^4;
        D4y(1, 1) = 7 / dy^4;   D4y(my, my) = 7 / dy^4;
    otherwise
        error('plate_matrices: unknown boundary condition "%s"', bc);
end

L = kron(Ix, D2y) + kron(D2x, Iy);
B = kron(Ix, D4y) + 2 * kron(D2x, D2y) + kron(D4x, Iy);
end