function pol = hermite( x, y, d, u)
	n = length(x)
	upperV = repmat(x', 1, 2*n)
	upperV = bsxfun(@power, upperV, [2*n-1:-1:0])
	lowerV = repmat(x', 1, 2*n)
	lowerV = bsxfun(@power, lowerV, [2*n-2:-1:0 , 1])
	lowerV = bsxfun(@times, lowerV, [2*n-1:-1:1 , 0])
	V = [upperV; lowerV]
	a = V \ [y.' ; d.' ]
	pol = polyval(a, u)
endfunction