%
%	Function yi = Interpol(x,y,xi);
%
function yi = interpol(x,y,xi);
%
[m,nx] = size(x);
[m,ni] = size(xi);
dx = x(2) - x(1);
%
q = floor((xi - x(1)) / dx) + 1;
qq = ceil((xi - x(1)) / dx);
p = max(q,2) + min(q,nx-1) - q;
s = max(0,sign(q)) + max(0,sign(nx - qq)) - 1;
%
b = [0,y(2:nx-1) - y(1:nx-2),0] / dx;
a = [0,y(3:nx) - 2 * y(2:nx-1) + y(1:nx-2),0] / (2 * dx^2);
%
yi = y(p-1) + (xi - x(p-1)) .* (b(p) + a(p) .* (xi - x(p)));
yi = s .* yi;
%
