% Sin curve using Cubic Spline Interpolation Method ...
% Author : Saket Soni
% Description : Numerial Assignment for plotting a sin curve using
%               using a Cubic Spline Interpolation Method.

x1 = 0;                 %   boundary values
xn = 15;                %

n = 15;                 %   no of actual values
en = 5;                 %   no of spline values to be generated between each pair of points

h = (xn-x1)/(n-1);

x = x1 : h : xn;        %   actual values x
y = sin(x);             %   actual values y


ex = zeros(1,(n-1)*en+n);       %   x for spline values
ey = zeros(1,(n-1)*en+n);       %   y for spline values

a = zeros(1,n-1);       %   for coefficients of degree zero
b = zeros(1,n-1);       %   for coefficients of degree one
c = zeros(1,n-1);       %   for coefficients of degree two
d = zeros(1,n-1);       %   for coefficients of degree three
cm = zeros(n-2);          %   matrix used for finding the C's.
cmr = zeros(1,n-2);       %   vector used for finding the C's

z = 0;                  %   for printing the x axis
z1 = 1;                 %   for printing the 1 line
zm1 = -1;               %   for printing the -1 line


% Now finding the C's ...
h = x(3) - x(2);
hl = x(2) - x(1);
cm(1,1) = 2*(h+hl);
cm(1,2) = h;
cmr(1) = 3*((y(3)-y(2))/h - (y(2)-y(1))/hl);
h = x(n) - x(n-1);
hl = x(n-1) - x(n-2);
cm(n-2,n-3) = hl;
cm(n-2,n-2) = 2*(h+hl);
cmr(n-2) = 3*((y(n)-y(n-1))/h - (y(n-1)-y(n-2))/hl);
for i = 3 : n-2
    h = x(i+1) - x(i);
    hl = x(i) - x(i-1);
    cm(i-1,i-2) = hl;
    cm(i-1,i-1) = 2*(h+hl);
    cm(i-1,i) = h;
    cmr(i-1) = 3*((y(i+1)-y(i))/h - (y(i)-y(i-1))/hl);
end

c(2 : n-1) = inv(cm)*transpose(cmr);
c(1) = 0;
c(n) = 0;


% Now finding the other coefficients ...
hl = x(2) - x(1);
a(1) = y(1);
d(1) = (c(2) - c(1))/(3*hl);
b(1) = (y(2) - y(1))/hl - d(1)*hl*hl;
for k = 2 : n-1
    a(k) = y(k);
    h = x(k+1) - x(k);
    d(k) = (c(k+1)-c(k))/(3*h);
    b(k) = b(k-1) + ( c(k) + c(k-1) ) * hl;
    hl = h;
end


% Now finding the spline values ...
j = 1;
for k = 1 : n-1
    ex(j) = x(k);
    ey(j) = y(k);
    j = j + 1;
    h = (x(k+1) - x(k))/(en+1);
    for i = 1 : en
        ex(j) = x(k) + h*i;
        t = ex(j) - x(k);
        t2 = t*t;
        ey(j) = a(k) + b(k)*t + c(k)*t2 + d(k)*t2*t;
        j = j + 1;
    end
end
ex((n-1)*en+n)=x(n);
ey((n-1)*en+n)=y(n);


% Now plotting the data ...
plot(x,y,ex,ey,ex,z,ex,z1,ex,zm1);

end