clear; clc;

% 输入参数：
% vec_x：已知数据点的x坐标（行向量）
% vec_y：已知数据点的y坐标（行向量，与vec_x长度相同）
% y01：左端点（vec_x(1)）的一阶导数
% yn1：右端点（vec_x(end)）的一阶导数
vec_x = [1, 2, 3, 4];  % 示例x坐标
vec_y = [1, 4, 9, 16]; % 示例y坐标
y01 = 2;               % 左端点一阶导数
yn1 = 8;               % 右端点一阶导数


[m, n] = size(vec_x);
h = zeros(1, n - 1);
h(1) = vec_x(2) - vec_x(1);
miu = zeros(1, n);
lambt = zeros(1, n);
MatrMomem = zeros(n);
d = zeros(1, n);

for j = 2 : n - 1
    h(j) = vec_x(j + 1) - vec_x(j);
    miu(j) = h(j - 1) / (h(j - 1) + h(j));
    lambt(j) = 1 - miu(j);
    d(j) = 6 * ((vec_y(j + 1) - vec_y(j)) / h(j) - (vec_y(j) - vec_y(j - 1)) / h(j - 1)) / (h(j - 1) + h(j));
    MatrMomem(j, j - 1 : j + 1) = [miu(j), 2, lambt(j)];
end

d(1) = 6 * ((vec_y(2) - vec_y(1)) / h(1) - y01) / h(1);
d(n) = 6 * (yn1 - (vec_y(n) - vec_y(n - 1)) / h(n - 1)) / h(n - 1);

MatrMomem(1, 1 : 2) = [2, 1];
MatrMomem(n, n - 1 : n) = [1, 2];
M = MatrMomem \ d';

% 输出结果：
% MatrMomem：三弯矩法的系数矩阵
% d：三弯矩法方程组的右端项向量
% M：求解得到的弯矩向量
MatrMomem, d, M