function [v,l,B] = solveEigenvalue(self,eigenvalueRange,coefstruct,tsecondOrder)
% solveEigenvalue - Internal function to solve eigenvalue problem
%

% Copyright 2021-2022 The MathWorks, Inc.

validateattributes(eigenvalueRange,{'numeric'},{'real', 'nonsparse',...
                                                'vector','numel',2,'nondecreasing'},2);

% Compute BC matrices, this would account for zero displacements.
[p,e,t] = self.Mesh.meshToPet();
u0 = pdeuxpd(p,0,self.PDESystemSize);
tdummy = 0:1;

femat = pde.DynamicDiscretizedPDEModel(self,p,e,t,coefstruct,u0,tdummy,tsecondOrder);

if (femat.vc || femat.vd || femat.va || femat.vq || femat.vh)
    error(message('pde:pdeModel:coefsNonlinear'))
end

K = femat.K;
A = femat.A;
Q = femat.Q;
H = femat.H;
M = femat.Mass;

% Impose Dirichlet BC.
[B,~]=pdenullorth(H);
K = B'*(K+Q+A)*B;
M = B'*M*B;

% If the asymmetry in the modal system is on the order of 10 times the
% numerical round-off error, then make the matrices symmetric.
K = self.checkForSymmetry(K);
M = self.checkForSymmetry(M);

opts.FrequencyRange = [eigenvalueRange(1),eigenvalueRange(2)];

if ~isempty(self.SolverOptions.MaxShift)
    opts.MaxShift = self.SolverOptions.MaxShift;
end

if ~isempty(self.SolverOptions.BlockSize)
    opts.BlockSize = self.SolverOptions.BlockSize;
end

[v,l] = matlab.internal.math.lanczos(K,M,opts);

%Recover Dirichlet BCs
if ~isempty(v)
    v=B*v;
end

end