R=0;
while R < 3
    prompt = 'Number of points? ';
    R=input(prompt);
end

x=[];
y=[];

[X,Y]=ginput(R);  
%%X=floor(X);
%%Y=floor(Y);
x = [X];
y = [Y];
n = length(x);

disp(x);
disp(y);



for i=1:n-1
    h(i) = x(i+1) - x(i);
end

for i=2:n-1
    al(i)=((3/h(i))*(y(i+1)-y(i)))-((3/h(i-1))*(y(i)-y(i-1)));
end

l(1)=1;

m(1)=0;

z(1)=0;

for i=2:n-1
    l(i)=(2*(x(i+1)-x(i-1)))-(h(i-1)*m(i-1));
    m(i)=h(i)/l(i);
    z(i)=(al(i)-h(i-1)*z(i-1))/l(i);
end

l(n)=1;

z(n)=0;

c(n)=0;

for j=n-1:-1:1  
    c(j)=z(j)-(m(j)*c(j+1));
    b(j)=(y(j+1)-y(j))/h(j)-(h(j)*(c(j+1)+2*c(j)))/3;
    d(j)=(c(j+1)-c(j))/(3*h(j));
end

for i=1:n-1
    a(i)=y(i);
end

p=x;

for k=1:n-1
    fun=@(x) a(k)+(b(k)*(x-p(k)))+(c(k)*(x-p(k)).^2)+(d(k)*(x-p(k)).^3) ;
    fplot(fun,[p(k),p(k+1)])
    hold on
end

axis([0,1,0,1]);