function [S] = hw8afn(m,x,y,z)

   h = x(2:m)-x(1:m-1);
   A = diag(2*(h(2:m-1)+h(1:m-2)))+diag(h(2:m-2),1)+diag(h(2:m-2),-1);
   b = (6*(y(2:m)-y(1:m-1))./h(1:m-1))';
   v = b(2:end)-b(1:end-1);
   zcoef(2:m-1) = (A\v)';
   zcoef(1) = 0;
   zcoef(m) = 0;

   i = 1;
   while x(i) < z
      i = i+1;
   end
   i = i-1;
   if i < 1
      i = 1;
   end
   S = zcoef(i)/(6*h(i))*(x(i+1)-z)^3+...
       zcoef(i+1)/(6*h(i))*(z-x(i))^3+...
       (y(i+1)/h(i)-zcoef(i+1)*h(i)/6)*(z-x(i))+...
       (y(i)/h(i)-zcoef(i)*h(i)/6)*(x(i+1)-z);

end
