Octave / square_interp

function y1 = square_interp (x, y, new_x)
y1 = [];
for i = 1 : length(new_x)
k = 3;
while (new_x(i) > x(k))
k += 1;
endwhile
a = (y(k) - y(k-2) - (y(k-1) - y(k-2))/(x(k-1) - x(k-2))*(x(k) - x(k-2))) / (x(k)^2 - (x(k-1) + x(k-2)) * x(k) - x(k-2)^2 + (x(k-1) + x(k-2)) * x(k-2));
b = (y(k-1) - y(k-2)) / (x(k-1) - x(k-2)) - a * (x(k-1) + x(k-2));
c = y(k-2) - a * x(k-2)^2 - b * x(k-2);
res = @(x) a * x .^2 + b * x + c;
y1 = [y1, res(new_x(i))];
endfor
endfunction