% plate_c.m -- Eigenvalue problem of rectangular plates with clamped edges
%  written by Hirano Yuuki, Meiji University, Feb 2012
%
%  \triangle^2 u=\lambda u
%  u=0, \partial u/\partial n=0
%  0<=x<=a, 0<=y<=b
%
% usage:
%   A=plate_c(1,1,640,640)
%   [v,d]=eigs(A,200,0)
%   plot_d(v(:,201-1),640,640)

function A=plate_c(a,b,nx,ny)
        hx=a/nx;
        hy=b/ny;
        m=nx-1;
        n=ny-1;
        Im=speye(m,m);
        In=speye(n,n);
        Jm=sparse(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1));
        Jn=sparse(diag(ones(n-1,1),1)+diag(ones(n-1,1),-1));
        Pm=sparse(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2));
        Pn=sparse(diag(ones(n-2,1),2)+diag(ones(n-2,1),-2));
        Am=-6*Im+4*Jm-Pm;
        Am(1,1)=-7;
        Am(m,m)=-7;
        An=-6*In+4*Jn-Pn;
        An(1,1)=-7;
        An(n,n)=-7;
        Bm=2*Im-Jm;
        Bn=2*In-Jn;

        A=-(kron(In,Am)/(hx*hx*hx*hx)+2*kron(-Bn,Bm)/(hx*hx*hy*hy)+kron(An,Im)/(hy*hy*hy*hy));

end
