Paste code

function NNatural = NNatural(X, Y)
    n = size(X, 2);
    for i = 1:(n-1)
        H(i) = X(i+1) - X(i);
        D(i) = (Y(i+1) - Y(i)) / H(i);
    end
    for i = 2:(n-1)
        U(i-1) = 6 * (D(i) - D(i-1));
    end
    
    A = zeros(n, n);
    for i = 1:(n-2)
        for k = 1:(i-1)
            A(i+1, k) = 0;
        end        
        A(i+1, i) = H(i);
        A(i+1, i+1) = 2 * (H(i) + H(i+1));
        A(i+1, i+2) = H(i+1);
        for k = (i+3):n
            A(i+1, k) = 0;
        end
    end
    A(:,n) = [];
    A(:,1) = [];
    A(n,:) = [];
    A(1,:) = [];
    
    M = SolveLUMatrix(A, U');
    M = horzcat([0], M')';
    M = horzcat(M', [0])';
    
    syms x;    
    hold on;
    for i = 1:(n-1)
        w = x - X(i);
        S(1) = Y(i);
        S(2) = D(i) - (H(i) * (2 * M(i) + M(i+1))) / 6;
        S(3) = M(i) / 2;
        S(4) = (M(i+1) - M(i)) / (6 * H(i));
        S = ((S(4) * w + S(3)) * w + S(2)) * w + Y(i);
        ezplot(S, [X(i), X(i+1)]);
    end
    axis([X(1) X(n) 0 10]);
end