Reconciling solvepde and assembleFEMatrices in PDE toolbox

I am trying to follow Matlab's example on dynamic analysis of a clamped beam . My objective is to extract the finite element matrices using assembleFEMatrices for this example, obtain the solution by time-stepping such as backward Euler, and comparing the result with Matlab's solvepde.

Let "model" be the model in the Matlab documentation. I extracted the FE matrices using:

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

Following the documentation on FE method and writing the PDE, if are the displacements in the directions, the resulting system of ODEs must be where z = [x ; y]. I can then rewrite this second order system into a first order as shown here .

However, I am unable to reconcile my solution with what I get from Matlab's solvepde.

My code is attached below. I changed some constants from the Matlab example so that the code run faster. I also considered the full beam instead of just half of it. I have also tried using lsim to time step the resulting system but it is somewhat comparable to what I get from backward euler.

Comments/suggestions much appreciated.

% 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