第一个猜测是使用前向欧拉方法。第一个雅可比使用有限差分。然后使用 NR 方法求解下一次迭代，并使用 Broyden 方法更新雅可比，直到满足容差。该图显示了我得到的结果与正确结果的比较。我可以对我的代码进行哪些改进，以便我朝着正确的方向前进（至少）。这是MATLAB代码。

function [t,y] = bdf1_qn_broyden(fcn,tspan,Y0,h,tol)

t = (tspan(1):h:tspan(2))';     %timestep vector
n = length(t);                  %no of steps
y = zeros(n,length(Y0));        %creting output array
y(1,:) = Y0;                    %first row of o/p is initial condition
i = 2;                              
jac = zeros(length(Y0));
unity = eye(length(Y0));

while i <= n 
    H = h*unity;
    for j = 1:1:length(Y0)
        jac(:,j) = fcn(t(i),y(i-1,:)' + H(:,j)) - fcn(t(i),y(i-1,:)');
    end
    iterate = (unity - jac)\unity;
    jac = (1/h)*jac;
    ycurr = y(i-1,:)' + h*fcn(t(i-1),y(i-1,:)');
    count = 0;
    diff = 1;
    while any(diff > tol) && (count < 25)
        curreval = fcn(t(i),ycurr);
        ynext = ycurr - iterate*(ycurr - y(i-1,:)' - h*curreval);
        delta = ynext - ycurr;
        deltaeval = fcn(t(i),ynext) - curreval;
        jac = jac + (1/dot(delta,delta))*(deltaeval - jac*delta)*delta';
        iterate = (unity - h*jac)\unity;
        diff = abs(delta);
        ycurr = ynext;
        count = count +1;
    end
    if count >= 25
    disp('Iterative method failed to converge within 25 steps');
    end
    y(i,:) = ynext';
    i = i+1;    
end
end

正如您可能已经猜到的那样，我是新手，我的代码效率很低。但是，如果我要慢一点，我希望至少是准确和慢的。如果不使用自适应步长，这是不可能的吗？BDF1 是否无法解决此类问题？我应该选择不同的方法吗？还是我应该更好地解决我的隐式函数？我的实施中有错误吗？我不知道该走哪条路。请指出我正确的方向。非常感谢！