clear; clc; close all;

% =========================
% 模版点：stencil
% =========================
k = -3:3; 
n = length(k);

% =========================
% 导数阶数
% =========================
m = 2;   

% =========================
% 构造线性系统（关键）
% =========================
A = sym(zeros(n));

for row = 1:n
    order = row - 1;
    A(row,:) = sym(k).^order / factorial(order);
end

b = sym(zeros(n,1));
b(m+1) = sym(1);   % 👈 第 m 阶导数

a = simplify(A \ b);

disp('符号解：');
disp(a);

% =========================
% 以下仅用于显示，提取统一分母
% =========================

den_list = [];

for i = 1:n
    if a(i) ~= 0
        [~,D] = numden(a(i));
        den_list = [den_list, double(D)];
    end
end

if isempty(den_list)
    common_den = sym(1);
else
    common_den = sym(den_list(1));
    for i = 2:length(den_list)
        common_den = lcm(common_den, sym(den_list(i)));
    end
end

a_int = simplify(a * common_den);

% =========================
% 拼分子
% =========================

num_str = "";

for i = 1:n
    if a_int(i) ~= 0
        coeff = a_int(i);
        shift = k(i);

        % 系数
        if coeff == 1
            c_str = "";
        elseif coeff == -1
            c_str = "-";
        else
            c_str = string(coeff);
        end

        % 下标
        if shift == 0
            u_str = "u_i";
        elseif shift > 0
            u_str = "u_{i+" + string(shift) + "}";
        else
            u_str = "u_{i" + string(shift) + "}";
        end

        term = c_str + u_str;

        if strlength(num_str) == 0
            num_str = term;
        else
            if startsWith(term, "-")
                num_str = num_str + term;
            else
                num_str = num_str + "+" + term;
            end
        end
    end
end

% =========================
% 分母处理
% =========================

if common_den == 1
    den_str = "(\Delta x)^" + string(m);
else
    den_str = string(common_den) + "(\Delta x)^" + string(m);
end

% =========================
% 自动导数符号
% =========================

if m == 1
    deriv_str = "\frac{\partial u}{\partial x}";
elseif m == 2
    deriv_str = "\frac{\partial^2 u}{\partial x^2}";
else
    deriv_str = "\frac{\partial^{" + string(m) + "} u}{\partial x^{" + string(m) + "}}";
end

% =========================
%  自动判断截断误差阶
% =========================

order = [];
for p = m+1 : m+10   % 往后多检查几阶
    moment = simplify(sum(a .* (sym(k)').^p / factorial(p)));
    if moment ~= 0
        order = p - m;
        break;
    end
end

% =========================
% 最终公式
% =========================

final_str = "\left(" + deriv_str + "\right)_i = \frac{" ...
    + num_str + "}{" + den_str + "} + O(\Delta x^" + string(order) + ")";

disp(' ');
disp(final_str);

% =========================
% 展示
% =========================
figure('Color','k');
text(0.05,0.5,['$' char(final_str) '$'], ...
    'Interpreter','latex', ...
    'Color','w', ...
    'FontSize',22);
axis off;