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: spm_ness_Lap.m
    path of file: ./repos/spm
    
    the code of the file until where you have to start completion: function [NESS] = spm_ness_Lap(M,x)
% Nonequilibrium steady-state under a Helmholtz decomposition
% FORMAT [NESS] = spm_ness_Lap(M,x)
%--------------------------------------------------------------------------
% M   - model specification structure
% Required fields:
%    M.f   - dx/dt   = f(x,u,P)  {function string or m-file}
%    M.pE  - P       = parameters of equation of motion
%    M.x   - (n x 1) = x(0) = expansion point
%    M.W   - (n x n) - precision matrix of random fluctuations
% x    - cell array of vectors specifying evaluation grid
%
% p0      - nonequilibrium steady-state
% X       - evaluation points of state space
% F       - expected flow
% f       - original flow
%
% NESS.H  - expected Hessian
% NESS.J  - expected Jacobian
% NESS.E  - Lyapunov exponents
% NESS.H2 - expected Euclidean norm of Hessian
% NESS.J2 - expected Euclidean norm of Jacobian
% NESS.D2 - correlation dimension
% NESS.bS - p0 = spm_softmax(spm_dctmtx(nx,nb)*bS);
% NESS.nb - number of basis functions
%__________________________________________________________________________

% Karl Friston
% Copyright (C) 2008-2022 Wellcome Centre for Human Neuroimaging


    
function [NESS] = spm_ness_Lap(M,x)
% Nonequilibrium steady-state under a Helmholtz decomposition
% FORMAT [NESS] = spm_ness_Lap(M,x)
%--------------------------------------------------------------------------
% M   - model specification structure
% Required fields:
%    M.f   - dx/dt   = f(x,u,P)  {function string or m-file}
%    M.pE  - P       = parameters of equation of motion
%    M.x   - (n x 1) = x(0) = expansion point
%    M.W   - (n x n) - precision matrix of random fluctuations
% x    - cell array of vectors specifying evaluation grid
%
% p0      - nonequilibrium steady-state
% X       - evaluation points of state space
% F       - expected flow
% f       - original flow
%
% NESS.H  - expected Hessian
% NESS.J  - expected Jacobian
% NESS.E  - Lyapunov exponents
% NESS.H2 - expected Euclidean norm of Hessian
% NESS.J2 - expected Euclidean norm of Jacobian
% NESS.D2 - correlation dimension
% NESS.bS - p0 = spm_softmax(spm_dctmtx(nx,nb)*bS);
% NESS.nb - number of basis functions
%__________________________________________________________________________

% Karl Friston
% Copyright (C) 2008-2022 Wellcome Centre for Human Neuroimaging


% get domain of phase-space and polynomial basis set
%==========================================================================
U      = spm_ness_U(M,x);      % get state space and flow
f      = U.f';                 % target flow
[nX,n] = size(U.X);            % number of sample points and states 

% get constraints (allowing G to be free parameters)
%--------------------------------------------------------------------------
J      = any(U.J,3);
J      = J'.*J;
[ks,kq,kg] = spm_NESS_constraints(U.o,J,M.K,M.L);
k      = find(~(kq | kg));
k      = find(~(kg));
nb     = size(U.b,2);
nQ     = n*n/2 + n/2;

pE.Qp  = zeros(nb*nQ,1);       % polynomial coefficients for solenoidal operator
pE.Rp  = zeros(1 ,n);          % polynomial coefficients for mean
pE.Sp  = zeros(nb,n,n);        % polynomial coefficients for Kernel

pC.Qp  = pE.Qp;
pC.Rp  = pE.Rp;
pC.Sp  = pE.Sp;

% constraints on solenoidal operator
%--------------------------------------------------------------------------
pC.Qp(k) = 1;

% constraints on mean
%--------------------------------------------------------------------------
for i = 1:n
    pE.Rp(1,i)  = mean(x{i});
    pC.Rp(1,i)  = 1;
end