You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: nrtHmxVibroSlab2d.m
    path of file: ./repos/Swan/gypsilabModified/nonRegressionTest/vibroAcoustic
    
    the code of the file until where you have to start completion: %+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal, Marc Bakry (c) 2017-2019.                 |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published
    
function (50% full, 10% decrease) 
    cutoff = vibsCutoff(1,L/5,L/10);
    
    % Green kernel function
    Gxy         = @(X,Y) femGreenKernel(X,Y,'[H0(kr)]',k0(i));
    gradyGxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]1',k0(i));
    gradyGxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]2',k0(i));
    gradyGxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[H0(kr)]3',k0(i));

    % Plane wave function
    PW         = @(X) exp(1i*k0(i)*X*X0') .* cutoff(X);
    gradxPW{1} = @(X) 1i*k0(i)*X0(1) .* PW(X);
    gradxPW{2} = @(X) 1i*k0(i)*X0(2) .* PW(X);
    gradxPW{3} = @(X) 1i*k0(i)*X0(3) .* PW(X);
    
    % Coupling coeff for Brackage-Werner simulation  
    beta = 1i*k0(i);
    
    % Quadrature and finite elements (volumn)
    omega = dom(mesh,3);
    U     = fem(mesh,'P1');

    % Quadrature and finite elements (boundary)
    sigma = dom(bound,3);
    u     = fem(bound,'P1');
    
    % Left-hand side
    tic
    [A,B,C,D] = vibsHmxBlockOperator(omega,U,sigma,u,cL,cT,rhoS,c0,rho0,f(i),tol);
    toc
    
    % Add dirichlet condition to x unknows (penalization)
    A(sub2ind(size(A),1:length(U),1:length(U))) = 1e15;
    
    % Right-hand side
    V    = cell(3,1);
    V{1} = - integral(sigma,ntimes(U,1),PW);
    V{2} = - integral(sigma,ntimes(U,2),PW);
    V{3} = integral(sigma,ntimes(u),gradxPW);
    
    % Resolution with Schur complement
    Fa     = decomposition(A);
    LHS    = @(V) D*V - C*(Fa \ (B.Ml*(B.Mr*V)) );      
    RHS    = V{end