function [ ev_y ] = spline_interpol( ev_x, x, y )

L=size(x,2)-1;
h=x(1,2:L+1)-x(1,1:L);
A=1/6*(2*diag(h(1,1:L-1))+2*diag(h(1,2:L))+diag(h(1,2:L-1),-1)+diag(h(1,2:L-1),1));
%display(A)
d=(y(1,3:L+1)-y(1,2:L))./h(1,2:L)-(y(1,2:L)-y(1,1:L-1))./h(1,1:L-1);
%display(d)
gamma=linsolve(A,d');
gamma=[0,gamma',0];
%display(gamma)
yd=(y(2:L+1)-y(1:L))./h(1,1:L)+h(1,1:L)/6.*(2*gamma(1,2:L+1)+gamma(1,1:L));
temp=yd(1,1)-h(1,1)/2*(gamma(1,1)+gamma(1,2));
yd=[temp,yd];
N=size(ev_x,2);
ev_y=zeros(size(ev_x));
for k=1:N
    for i=2:L+1
        if ev_x(1,k)<x(1,i)
            ev_y(1,k)=y(1,i)+yd(1,i)*(ev_x(1,k)-x(1,i))+gamma(1,i)/2*(ev_x(1,k)-x(1,i))^2+(gamma(1,i)-gamma(1,i-1))/(6*h(1,i-1))*(ev_x(1,k)-x(1,i))^3;
            break;
        end
    end
end
