%es 2 19.09.2022
clear all
close all
clc

x=linspace(0,3,10)';
f= @(x) -6*sin(2*x) + x.^2 +7*x;
y=f(x);

xx=linspace(0,3);

h=diff(x);
v=diff(y)./h;

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

%Az=b
A=diag(D)+diag(C,-1)+diag(C,1);
b=6*(v(2:end)-v(1:end-1));
z=A\b;
z=[0; z; 0];

%Costruisco yy
yy=zeros(size(xx));

for ii=1:length(x)-1
    I=(x(ii)<=xx)&(x(ii+1)>=xx);
    
    C=y(ii+1)./h(ii)-h(ii)/6*z(ii+1);
    
    D=y(ii)/h(ii)-h(ii)/6*z(ii);
    yy=yy+(z(ii+1)/6/h(ii).*(xx-x(ii)).^3+z(ii)/6/h(ii).*(x(ii+1)-xx).^3+...
        +C*(xx-x(ii))+D*(x(ii+1)-xx)).*I;
end

figure
plot(x,y,'*')
hold on
plot(xx,yy)