function yy = cubicspline(xx,x,y)

n  = size(x       ,2);
nn = size(xx      ,2);
A  = zeros(4*(n-1)  ); 
b  = zeros(4*(n-1),1);
yy = zeros(nn     ,1);
for i = 1:n-1
  %            fits f_i
  A(4*(i-1)+1,4*(i-1)+  [ 1  2  3 4]) =  [  x(i  )^3   x(i  )^2 x(i  )^1 x(i  )^0];   b(4*(i-1)+1,1) = y(i  );
  % f          its f_i+1
  A(4*(i-1)+2,4*(i-1)+  [ 1  2  3 4]) =  [  x(i+1)^3   x(i+1)^2 x(i+1)^1 x(i+1)^0];   b(4*(i-1)+2,1) = y(i+1);
  if i == 1  % 2 derivative x_i+1 for last interval
    A(4*(i-1)+3,4*(n-1)+[ 1  2  3 4]) =  [6*x(i+1)^1 2*x(i+1)^0 0        0       ];   b(4*(i-1)+3,1) = 0     ;
  else       % 1 derivative x_i
    A(4*(i-1)+3,4*(i-1)+[ 1  2  3 4]) =  [3*x(i  )^2 2*x(i  )^1 x(i  )^0 0       ];   b(4*(i-1)+3,1) = 0;
    A(4*(i-1)+3,4*(i-1)+[-3 -2 -1 0]) = -[3*x(i  )^2 2*x(i  )^1 x(i  )^0 0       ];
  end
  %            2 derivative x_i
  if i == n-1
    A(4*(i-1)+4,4*(i-1)+[ 1  2  3 4]) =  [6*x(i  )^1 2*x(i  )^0 0        0       ];   b(4*(i-1)+4,1) = 0     ;
  else
    A(4*(i-1)+4,4*(i-1)+[ 1  2  3 4]) =  [3*x(i  )^2 2*x(i  )^1 x(i  )^0 0       ];   b(4*(i-1)+4,1) = 0;
    A(4*(i-1)+4,4*(i-1)+[ 5  6  7 8]) = -[3*x(i  )^2 2*x(i  )^1 x(i  )^0 0       ];
  end
end
al = A\b;
icnt = 1;
for ii = 1:nn
  if xx(ii) > x(min(n-1,icnt + 1)),     icnt = min(n-1,icnt + 1);   end
  abcd = al(4*(icnt-1)+[1 2 3 4],1);
  yy(ii) = abcd(1)*xx(ii)^3 + abcd(2)*xx(ii)^2 + abcd(3)*xx(ii)^1 + abcd(4)*xx(ii)^0;
end
