 function [M,K] = beam_mats(N)

% function [M,K] = beam_mats(N)
%
% constructs the mass and stiffness matrices associated with the operator
% Lu = u'''' for x in [0,1] with u(0)=u''(0)=u(1)=u''(1)=0
% using piecewise cubic Hermite basis functions.  [Burkhardt]

 h = 1/(N+1);

% build the mass matrix, M

 A = h*[2/105 13/420; 13/420 26/35];
 B = h*[-1/140 0; 0 9/70];
 C = h*[0 -13/420; 0 0]; 

 T1 = diag(ones(N,1),1)+diag(ones(N,1),-1);
 T2 = diag(ones(N-1,1),-2);
 M = kron(eye(N+1),A) + kron(T1,B) + kron(T2,C) + kron(T2',C');

% set the boundary conditions for M
 M(1,1) = h/105;
 M(2*N:2*N+1,2*N+2) = h*[-13/420; -1/140];
 M(2*N+2,[2*N:2*N+2]) = h*[-13/420 -1/140 1/105]; 

% build the stiffness matrix, K

 A = [8 -6; -6 24]/(h^3);
 B = [2 0; 0 -12]/(h^3);
 C = [0 6; 0 0]/(h^3); 

 T1 = diag(ones(N,1),1)+diag(ones(N,1),-1);
 T2 = diag(ones(N-1,1),-2);
 K = kron(eye(N+1),A) + kron(T1,B) + kron(T2,C) + kron(T2',C');
 
% set the boundary conditions for K
 K(1,1) = 4/(h^3);
 K(2*N:2*N+1,2*N+2) =   [6; 2]/(h^3);
 K(2*N+2,[2*N:2*N+2]) = [6 2 4]/(h^3);
