%Spline naturali
% derivate seconde negli estremi nulle;

clear all
close all
clc

x=[];
y=[];
xx=linspace(x(1),x(end));
[yy] = SPLINENAT(x,y,xx);

f= @(xx) SPLINENAT(x,y,xx);
fplot(f,[x(1),x(end)])

function [yy]=SPLINENAT(x,y,xx)
h=diff(x);
v=diff(y)./h;

D=2*(x(3:end)-x(1:end-2));
C=h(2:end-1);              

A=diag(D)+diag(C,-1)+diag(C,1);
b=6*diff(v);
M=A\b;
M=[0; M; 0];

yy=zeros(size(xx));
for ii=1:length(x)-1

    if ii==1
        I=(x(ii)<=xx)&(x(ii+1)>=xx);
    else
        I=(x(ii)<xx)&(x(ii+1)>=xx);
    end
    
    A=M(ii+1)/6/h(ii).*(xx-x(ii)).^3;
    B=M(ii)/6/h(ii).*(x(ii+1)-xx).^3;
    C=y(ii+1)./h(ii)-h(ii)/6*M(ii+1);
    D=y(ii)/h(ii)-h(ii)/6*M(ii);
    
    yy=yy+(A+B+C*(xx-x(ii))+D*(x(ii+1)-xx)).*I;
end

end