Why am I "Out of memory" because of an "Error using lu"

I Need to solve this PDE (the way they recommend here: Solve PDE ), but I got the following Problem:

Here´s the full Code:

%%Model parameters

% https://de.mathworks.com/help/pde/ug/c-coefficient-for-systems-for-specifycoefficients.html#bu5tna0-1

global N

N = 3;

global nu

nu = 0.34;

global E

E = 70;

global rho

rho = 2.6989;

global A

A = rho * (1- 2 * nu) * (1 + nu) / E;

global B

B = 1 - nu;

global C

   C = 1/2 - nu;

global D

D = 3/2 - 2 * nu;

% time grid

tlist = 0:0.025:0.1;

model = createpde(N);

% Import geometry into the container.

importGeometry(model,'zylinder.stl');

% View the geometry with face labels.

pdegplot(model,'FaceLabels','off')

%%Specify PDE Coefficients

% Include the PDE COefficients in |model|.

specifyCoefficients(model, 'm', @mCoeffFunction, 'd', 0, 'c', @cCoeffFunction, 'a', @aCoeffFunction, 'f', @fCoeffFunction);

% set initial conditions

u0 = [0.01; 0; 0];

ut0 = [0.1; 0; 0];

setInitialConditions(model,u0, ut0);

% Create zero Dirichlet boundary conditions on all faces.

applyBoundaryCondition(model,'dirichlet','face',1:3,'u',[0.01,0,0]); % noch falsche Randbedingung, da unklar wie einzugeben (Rb ist PDE).
% Create a mesh

generateMesh(model);

pdemesh(model)

solve the PDE

result = solvepde(model, tlist);

u = result.NodalSolution;
pdeplot3D(model,'ColorMapData',u(:,1))

title('u(1) --> u')

function mMatrix = mCoeffFunction(region, state)

      % https://de.mathworks.com/help/pde/ug/m-d-or-a-coefficient-for-systems.html#bu5xcw2-1

      % x is r here

      global N A B C D

      Nr = length(region.x);

      mMatrix = zeros(N^2, Nr);

      mMatrix(1, :) = -region.x .* A;

      mMatrix(5, :) = -region.x .* A;

      mMatrix(9, :) = -region.x .* A;

end

function dMatrix = dCoeffFunction(region, state)

      % https://de.mathworks.com/help/pde/ug/m-d-or-a-coefficient-for-systems.html#bu5xcw2-1

      global N A B C D

      Nr = length(region.x);

      dMatrix = zeros(N^2, Nr);

end

function cMatrix = cCoeffFunction(region, state)

      % https://de.mathworks.com/help/pde/ug/m-d-or-a-coefficient-for-systems.html#bu5xcw2-1

      global N A B C D

      Nr = length(region.x);

      cMatrix = zeros(9*N^2, Nr);

      cMatrix(getRowOfCMatrix(1, 1, 1, 1), :) = B .* region.x;

      cMatrix(getRowOfCMatrix(1, 1, 2, 2), :) = C ./ region.x;

      cMatrix(getRowOfCMatrix(1, 1, 3, 3), :) = C .* region.x;

      cMatrix(getRowOfCMatrix(1, 2, 1, 2), :) = 1/4;

      cMatrix(getRowOfCMatrix(1, 2, 2, 1), :) = 1/4;

      cMatrix(getRowOfCMatrix(1, 3, 1, 3), :) = 1/4 .* region.x;

      cMatrix(getRowOfCMatrix(1, 3, 3, 1), :) = 1/4 .* region.x;

      cMatrix(getRowOfCMatrix(2, 1, 1, 2), :) = 1/4;

      cMatrix(getRowOfCMatrix(2, 1, 2, 1), :) = 1/4;

      cMatrix(getRowOfCMatrix(2, 3, 3, 2), :) = 1/4;

      cMatrix(getRowOfCMatrix(2, 3, 2, 3), :) = 1/4;

      cMatrix(getRowOfCMatrix(2, 2, 1, 1), :) = C .* region.x;

      cMatrix(getRowOfCMatrix(2, 2, 2, 2), :) = B ./ region.x;

      cMatrix(getRowOfCMatrix(2, 2, 3, 3), :) = C .* region.x;

      cMatrix(getRowOfCMatrix(3, 1, 3, 1), :) = 1/4 .* region.x;

      cMatrix(getRowOfCMatrix(3, 1, 1, 3), :) = 1/4 .* region.x;

      cMatrix(getRowOfCMatrix(3, 3, 1, 1), :) = C .* region.x;

      cMatrix(getRowOfCMatrix(3, 3, 2, 2), :) = C ./ region.x;

      cMatrix(getRowOfCMatrix(3, 3, 3, 3), :) = B .* region.x;

end

function aMatrix = aCoeffFunction(region, state)

      % https://de.mathworks.com/help/pde/ug/m-d-or-a-coefficient-for-systems.html#bu5xcw2-1

      global N A B C D

      Nr = length(region.x);

      aMatrix = zeros(N^2, Nr);

      aMatrix(1, :) = -B ./ region.x;

      aMatrix(5, :) = -C ./ region.x;   

end

function fMatrix = fCoeffFunction(region, state)

      % https://de.mathworks.com/help/pde/ug/m-d-or-a-coefficient-for-systems.html#bu5xcw2-1

      % state.ux(1,:) --> ur, state.uy(2, :) --> vphi, ...

      global N A B C D

      Nr = length(region.x);

      fMatrix = zeros(N, Nr);

      fMatrix(1, :) = -B .* state.ux(1, :) + D * 1 ./ region.x .* state.uy(2, :);

      fMatrix(2, :) = -D ./ region.x .* state.uy(1, :) - C .* state.ux(2, :);

      fMatrix(3, :) = -1/2 .* state.uz(1, :) - C .* state.ux(3, :);

end

function row = getRowOfCMatrix(i, j, k, l)

      global N

      row = 9 * N * (j - 1) + 9 * i + 3 * l + k - 12;

end