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: nrtHmxMaxwellCFIE.m
    path of file: ./repos/Swan/gypsilabModified/nonRegressionTest/scattering3d
    
    the code of the file until where you have to start completion: %+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal & Francois
    
function --> G(x,y) = exp(ik|x-y|)/|x-y| 
Gxy    = @(X,Y) femGreenKernel(X,Y,'[exp(ikr)/r]',k);
Hxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]1',k) ;
Hxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]2',k) ;
Hxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]3',k) ;

% Finite elements
u = fem(sphere,'RWG');
v = fem(sphere,'RWG');

% Finite element mass matrix --> \int_Sx psi(x)' psi(x) dx
Id = integral(sigma,u,v);

% Finite element boundary operator
tic
T = 1i*k/(4*pi)*integral(sigma, sigma, v, Gxy, u, tol) ...
    - 1i/(4*pi*k)*integral(sigma, sigma, div(v), Gxy, div(u), tol) ;
toc

% Regularization
tic
Tr = 1i*k/(4*pi)*regularize(sigma, sigma, v, '[1/r]', u) ...
      - 1i/(4*pi*k)*regularize(sigma, sigma, div(v), '[1/r]', div(u));
T  = T + Tr; 
toc

% Finite element boundary operator
tic
nxK = 1/(4*pi) * integral(sigma, sigma, nx(v), Hxy, u, tol); 
toc

% Regularization
tic
nxKr = 1/(4*pi) * regularize(sigma, sigma, nx(v), 'grady[1/r]', u);
nxK  = nxK + nxKr;
toc

% Left hand side
LHS  = - beta * T  + (1-beta) * (0.5*Id - nxK);
LHSr = - beta * Tr + (1-beta) * (0.5*Id - nxKr);

% Right hand side
RHS = beta*integral(sigma,v,PWE) - (1-beta)*integral(sigma,nx(v),PWH);


%%% SOLVE LINEAR PROBLEM
disp('~~~~~~~~~~~~~ SOLVE LINEAR PROBLEM ~~~~~~~~~~~~~')

% ILU preconditionner
tic
[L,U] = ilu(LHSr);
toc

% Iterative solver
tic
J = gmres(@(V) LHS*V,RHS,[],tol,1000,L,U);
toc


%%% INFINITE SOLUTION
disp('~~~~~~~~~~~~~ INFINITE RADIATION ~~~~~~~~~~~~~')

% Plane waves direction
Ninf  = 1e3;
theta = 2*pi/1e3 .* (1:Ninf)';
nu    = [sin(theta),zeros(size(theta)),cos(theta)];

% Green kernel function
xdoty = @(X,Y) X(:,1).*Y(:,1) + X(:,2).*Y(:,2) + X(:,3).*Y(:,3); 
Ginf  = @(X,Y) exp(-1i*k*xdoty(X,Y));

% Finite element infinite operator --> \int_Sy exp(ik*nu.y) * psi(y) dx
Tinf = integral(nu,sigma,Ginf,v);
sol  = 1i*k/(4*pi)*cross(nu, cross([Tinf{1}*J, Tinf{2}*J, Tinf{3}*J], nu));

% Radiation infinie de reference, convention e^(+ikr)/r
nMax = 100; refInf = zeros(Ninf,1);
if E(1) == 1
    for jj = 1:Ninf
        refInf(jj,:) = sphereMaxwell(1, -f, theta(jj), 0.0, nMax);
    end