ClampedCubicSpline

clc
clear
x0=input('x0=');
y0=input('y0=');
FPO=input('FPO=');
FPN=input('FPN=');
x=input('x=');
n=length(x0);
for i=1:n
a(i)=y0(i);
end
for i=1:n-1
h(i)=x0(i+1)-x0(i);
end
g(1)=3*(a(2)-a(1))/h(1)-3*FPO;
g(n)=3*FPN-3*(a(n)-a(n-1))/h(n-1);
for i=2:n-1
g(i)=3*(a(i+1)-a(i))/(h(i))-3*(a(i)-a(i-1))/(h(i-1));
end
l(1)=2*h(1);
u(1)=0.5;
z(1)=g(1)/l(1);
for i=2:n-1
l(i)=2*(x(i+1)-x(i-1))-h(i-1)*u(i-1);
u(i)=h(i)/l(i);
z(i)=(g(i)-h(i-1)*z(i-1))/l(i);
end
l(n)=h(n-1)*(2-u(n-1));
z(n)=(g(n)-h(n-1)*z(n-1))/l(n);
c(n)=z(n);
for j=n-1:-1:1
c(j)=z(j)-u(j)*c(j+1);
b(j)=(a(j+1)-a(j))/h(j)-h(j)*(c(j+1)+2*c(j))/3;
d(j)=(c(j+1)-c(j))/(3*h(j));
end
for j=1:n-1
s(j)=a(j)+b(j)*(x-x0(j))+c(j)*(x-x0(j))^2+d(j)*(x-x0(j))^3;
end
disp(s)