Программа кубического сплайна

Я пытаюсь написать программу интерполяции кубическим сплайном. Я написал программу, но график не отображается правильно. Сплайн использует естественную границу условия (второй dervative в начальном / конечном узле равен 0). Код находится в Matlab и показан ниже,

clear all
%Function to Interpolate
k = 10;                    %Number of Support Nodes-1
xs(1) = -1;
for j = 1:k
    xs(j+1) = -1 +2*j/k;   %Support Nodes(Equidistant)
end;
fs = 1./(25.*xs.^2+1);     %Support Ordinates
x = [-0.99:2/(2*k):0.99];  %Places to Evaluate Function
fx = 1./(25.*x.^2+1);      %Function Evaluated at x

%Cubic Spline Code(Coefficients to Calculate 2nd Derivatives)

f(1) = 2*(xs(3)-xs(1));
g(1) = xs(3)-xs(2);
r(1) = (6/(xs(3)-xs(2)))*(fs(3)-fs(2)) + (6/(xs(2)-xs(1)))*(fs(1)-fs(2));
e(1) = 0;

for i = 2:k-2
    e(i) = xs(i+1)-xs(i);
    f(i) = 2*(xs(i+2)-xs(i));
    g(i) = xs(i+2)-xs(i+1);
    r(i) = (6/(xs(i+2)-xs(i+1)))*(fs(i+2)-fs(i+1)) + ...
           (6/(xs(i+1)-xs(i)))*(fs(i)-fs(i+1));
end
e(k-1) = xs(k)-xs(k-1);
f(k-1) = 2*(xs(k+1)-xs(k-1));
r(k-1) = (6/(xs(k+1)-xs(k)))*(fs(k+1)-fs(k)) + ...
         (6/(xs(k)-xs(k-1)))*(fs(k-1)-fs(k));

%Tridiagonal System

i = 1;
A = zeros(k-1,k-1);
while i < size(A)+1;
    A(i,i) = f(i);
    if i < size(A);
        A(i,i+1) = g(i);
        A(i+1,i) = e(i);
    end
    i = i+1;
end

for i = 2:k-1                             %Decomposition
    e(i) = e(i)/f(i-1);
    f(i) = f(i)-e(i)*g(i-1);
end

for i = 2:k-1                             %Forward Substitution 
    r(i) = r(i)-e(i)*r(i-1);
end

xn(k-1)= r(k-1)/f(k-1);
for i = k-2:-1:1                          %Back Substitution
    xn(i) = (r(i)-g(i)*xn(i+1))/f(i);
end

%Interpolation

if (max(xs) <= max(x))
    error('Outside Range'); 
end

if (min(xs) >= min(x))
    error('Outside Range'); 
end


P = zeros(size(length(x),length(x)));
i = 1;
for Counter = 1:length(x)
    for j = 1:k-1
        a(j) = x(Counter)- xs(j);
    end
    i = find(a == min(a(a>=0)));
    if i == 1
        c1 = 0;
        c2 = xn(1)/6/(xs(2)-xs(1));
        c3 = fs(1)/(xs(2)-xs(1));
        c4 = fs(2)/(xs(2)-xs(1))-xn(1)*(xs(2)-xs(1))/6;
        t1 = c1*(xs(2)-x(Counter))^3;
        t2 = c2*(x(Counter)-xs(1))^3;
        t3 = c3*(xs(2)-x(Counter));
        t4 = c4*(x(Counter)-xs(1));
        P(Counter) = t1 +t2 +t3 +t4;
    else
        if i < k-1
        c1 = xn(i-1+1)/6/(xs(i+1)-xs(i-1+1));
        c2 = xn(i+1)/6/(xs(i+1)-xs(i-1+1));
        c3 = fs(i-1+1)/(xs(i+1)-xs(i-1+1))-xn(i-1+1)*(xs(i+1)-xs(i-1+1))/6;
        c4 = fs(i+1)/(xs(i+1)-xs(i-1+1))-xn(i+1)*(xs(i+1)-xs(i-1+1))/6;
        t1 = c1*(xs(i+1)-x(Counter))^3;
        t2 = c2*(x(Counter)-xs(i-1+1))^3;
        t3 = c3*(xs(i+1)-x(Counter));
        t4 = c4*(x(Counter)-xs(i-1+1));
        P(Counter) = t1 +t2 +t3 +t4;
        else
        c1 = xn(i-1+1)/6/(xs(i+1)-xs(i-1+1));
        c2 = 0;
        c3 = fs(i-1+1)/(xs(i+1)-xs(i-1+1))-xn(i-1+1)*(xs(i+1)-xs(i-1+1))/6;
        c4 = fs(i+1)/(xs(i+1)-xs(i-1+1));
        t1 = c1*(xs(i+1)-x(Counter))^3;
        t2 = c2*(x(Counter)-xs(i-1+1))^3;
        t3 = c3*(xs(i+1)-x(Counter));
        t4 = c4*(x(Counter)-xs(i-1+1));
        P(Counter) = t1 +t2 +t3 +t4;
        end    
    end
end

P = P';
P(length(x)) = NaN;

plot(x,P,x,fx)

Когда я запускаю код, функция интерполяции не является симметричной и не сходится правильно. есть предложения о проблемах в моем коде? Спасибо.