function a = QuadSpline(x,y)

%

% x and y are equal length column vectors representing the points (x1,y1),...,(xn,yn)

%

% a is a vector [a11 a12 a13 a21 a22 a23 ... a(n-1)1 a(n-1)2 a(n-1)3]'

% such that:

%

% (1) qi(x) = ai1 + ai2(x-x_i) + ai3(x-x_i)(x-x_i+1) is the local cubic interpolant

% of (x_i,y_i) and (x_(i+1), y_(i+1)) on the interval [x_i, x_i+1] for i = 1,..,n-1

%

% (2) qi'(xi+1) = q_i+1'(xi+1)

%

% (3) The piecewise quadratic funtion has continuous second derivative at x_(n/2)

%

n = size(x,1);

V = zeros(3*n-3,3*n-3);

%Build the matrix for the linear system V*a = yBig

for i=1:n-1

V(3*(i-1)+1,3*(i-1)+1) = 1; %ai1 = yi

V(3*(i-1)+2,3*(i-1)+1) = 1; %ai1 + ai2*dx = y(i+1)

V(3*(i-1)+2,3*(i-1)+2) = x(i+1)-x(i); % ...for above eqn...

if(i<n-1)

%Encode constraint (2) n-2 times (once for each interior point xi).

V(3*(i-1)+3,3*(i-1)+2:3*(i-1)+6) = [1 (x(i+1)-x(i)) 0 -1 -(x(i+1)-x(i+2))];

else

%Encode constraint (3) once

V(3*n-3,3*ceil(n/2)-3) = 1;

V(3*n-3,3*ceil(n/2)) = -1;

end

end

%Build the right-hand-side vector yBig

yBig = zeros(3*n-3,1); %rhs for constraint (2) is zeros

yBig((1:3:3*(n-1)-2)') = y(1:n-1); %for ai1 equqtions

yBig((2:3:3*(n-1)-1)') = y(2:n); %for ai1 + ai2*dx equations

%Solve the linear system

a = V \ yBig;