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: lm.m
    path of file: ./repos/spinw/swfiles/+ndbase
    
    the code of the file until where you have to start completion: function [p0,yHat,stat] = lm(dat,func0,p0,varargin)
% least square refinement of parameters using Levenberg Marquardt method
%
% [pOpt,fVal,stat] = NDBASE.LM(dat,func,p0,'Option1','Value1',...)
%
% Levenberg Marquardt curve-fitting function, minimizes sum of weighted
% squared residuals.
%
% Input:
%
% dat       Data to be fitted stored in a structure with
    
function values.
% yDat  Vector of data to be fitted.
% weight The weighting vector for least squares fit ...
% dp 	Fractional increment of p for numerical derivatives:
%           dp(j)>0     central differences calculated
%           dp(j)<0     one sided differences calculated
%           dp(j)=0     sets corresponding partials to zero, i.e. holds p(j)
%                       fixed.
%       Default value is 1e-3.
%
% Output:
%
% JtWJ      Linearized Hessian matrix (inverse of covariance matrix).
% JtWdy   	Linearized fitting vector.
% chiSq     Chi-squared criteria, weighted sum of the squared residuals.
% yHat  	Model evaluated with parameters p.
% J         Jacobian matrix.
%

%   Henri Gavin, Dept. Civil & Environ. Engineering, Duke Univ. November 2005
%   modified from: ftp://fly.cnuce.cnr.it/pub/software/octave/leasqr/
%   Press, et al., Numerical Recipes, Cambridge Univ. Press, 1992, Chapter 15.

global   lm_iteration  lm_func_calls

% number of varying parameters
Np = numel(p);
% number of varying parameters
Nv = sum(dp~=0);

% evaluate model using parameters 'p'
yFunc = func(x,p');
lm_func_calls = lm_func_calls + 1;

if ~rem(lm_iteration,2*Np) || dX2 > 0
    % finite difference
    J = jacobian(func,x,p,yFunc,dp);
else
    % rank-1 update
    J = broyden(pOld,yOld,J,p,yFunc,dp);
end