Extracting FEM matrices in matlab pde toolbox

Question

I am trying to follow the dynamic linear elasticity in Matlab, link here. My goal is to extract the FE Matrices using the function assembleFEMatrices in matlab and solve the resulting system of second-order ODE's by Backward Euler, for example.

I am not having any success reconciling the result from solvepde with the result from my code.

Let ''model'' be the PDE model specified according to the Matlab example.

To extract the FE Matrices, I do:

FEM = assembleFEMatrices(model);
M = FEM.M;
K = FEM.K;
F = FEM.F;
G = FEM.G;

This is following Matlab's documentation on FEM, on writing the PDE, and on extracting the FE Matrices.

Let $x,y$ be the displacement in the $x,y$ directions, respectively. The resulting system of ODE's should be

$M \begin{bmatrix} \ddot{x} \\ \ddot{y} \end{bmatrix} + K \begin{bmatrix} x \\ y \end{bmatrix} = G+F$.

Here, $M$ is the mass matrix, $K$ is the stiffness matrix, $G$ is from the Neumann BC while $F$ is for the forcing term.

I then write the above system as

$ \begin{bmatrix} \ddot{x} \\ \ddot{y} \end{bmatrix} + M^{-1}K \begin{bmatrix} x \\ y \end{bmatrix} = M^{-1}(G+F)$.

Then to reduce this to a 1st-order ODE, I introduce $z_1 = \begin{bmatrix} x \\ y \end{bmatrix}$ and $z_2 = \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix}$

to obtain

$\frac{d}{dt} \begin{bmatrix} z_1 \\ z_2 \end{bmatrix} = \begin{bmatrix} 0 & I \\ -M^{-1}K & 0 \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \end{bmatrix} + \begin{bmatrix} 0 \\ M^{-1}(G+F) \end{bmatrix}$.

I then solve the above system using backward euler and compared $z_1$ with the result I get from Matlab's solvepde but they don't match. I don't think the problem is with my solution approach because I tried Matlab's lsim and I get almost the same results as backward euler. Did I incorrectly extract and misinterpret the FE Matrices?

Suggestions appreciated.

Edit: Code below. I slightly changed some constants in the Matlab example for the code to run faster. I also considered the full beam instead of just half the beam. In my code below, mySoln is the solution I get from assembleFEMatrices and transforming the 2nd order ODE to a first order ODE as I described above. results.NodalSolution is the solution obtained from Matlab's solvepde. I can't get this to reconcile with mySoln.

% Create PDE model

N = 2; % Two PDE in plane stress elasticity
model = createpde(N);
blength = 1.5; % Beam length, in.
height = .01; % Thickness of the beam, in.

l2 = blength/2;
h2 = height/2;
rect = [3 4 0 l2 l2 0 -h2 -h2  h2 h2]';
g = decsg(rect,'R1',('R1')');
pg = geometryFromEdges(model,g);

figure
pdegplot(g,'EdgeLabels','on');
title('Geometry With Edge Labels Displayed');
axis([-.1 1.1*l2 -5*h2 5*h2]);

E = 3.0e7; % Young's modulus of the material, lbs/in^2
gnu = .3; % Poisson's ratio of the material
rho = .3/386; % Density of the material
G = E/(2.*(1+gnu));
mu = 2*G*gnu/(1-gnu);
c = [2*G+mu; 0; G;   0; G; mu; 0;  G; 0; 2*G+mu];
f = [0 0]'; % No body forces
specifyCoefficients(model, 'm', rho, 'd', 0, 'c', c, 'a', 0, 'f', f);

% BCs
symBC = applyBoundaryCondition(model,'dirichlet','Edge',4,'u',[0 0]);
clampedBC = applyBoundaryCondition(model,'dirichlet','Edge',2,'u',[0 0]);
sigma = 5e1;
presBC = applyBoundaryCondition(model,'neumann','Edge',3,'g',[0 sigma]);
setInitialConditions(model,0,0);
generateMesh(model,'Hmax',height/2,'MesherVersion','R2013a');

% Solve pde using matlab pde toolbox
tfinal = 3e-3;
tlist = linspace(0,tfinal,100);
result = solvepde(model,tlist);

%%
% plot appearance of beam

figure;

X_nodes = result.Mesh.Nodes(1,:);
X_nodes = X_nodes(:);
Y_nodes = result.Mesh.Nodes(2,:);
Y_nodes = Y_nodes(:);

for k = 1:size(result.SolutionTimes,2)

    % Extract displacement at time k
    disp_X = result.NodalSolution(:,1,k);
    disp_X = disp_X(:);

    disp_Y = result.NodalSolution(:,2,k);
    disp_Y = disp_Y(:);

    plot(X_nodes + disp_X, Y_nodes + disp_Y,'*')
    pause(0.1)

end

%%

% try to recover solution by assembling the FEM matrix and time-stepping it
FEM = assembleFEMatrices(model);
Mmat = FEM.M;
Kmat = FEM.K;
Fmat = FEM.F;
Gmat = FEM.G;

nNodes = size(model.Mesh.Nodes,2);

Amat = [zeros(2*nNodes,2*nNodes) eye(2*nNodes); -Mmat\Kmat zeros(2*nNodes,2*nNodes)];
Bmat = [zeros(2*nNodes,1); Mmat\(Fmat + Gmat)];

% Do backward euler here

deltaT = tlist(2) - tlist(1);
IminushA = sparse(eye(4*nNodes)) - deltaT*Amat;

mySoln = zeros(4*nNodes,length(tlist));

for k = 2:length(tlist)
   fprintf('%d...\n',k)
   mySoln(:,k) = IminushA\(mySoln(:,k-1) + deltaT*Bmat);
end