unequalSimpson <- function(nodes, values)
{
	if(length(nodes) != length(values))
	{stop("Simpson: number of nodes unequal to number of values")}
	
	len <- length(nodes)
	pieces <- (len - 1)/2
	if(round(pieces) != pieces)
	{stop("Simpson: number of nodes must equal 2*n + 1")}


	#calculates simpsons rule for all pieces defined by three
	#successive points starting with the first one
	simpEvals1 <- vector(mode = "numeric", length = pieces)

	for(i in seq(1,pieces))
	{
		idx0 <- 2*i - 1
		idx1 <- 2*i
		idx2 <- 2*i + 1
		
		y0 <- values[idx0]
		y1 <- values[idx1]
		y2 <- values[idx2]
		
		delX1 <- nodes[idx1] - nodes[idx0]
		delX2 <- nodes[idx2] - nodes[idx1]	
	
		simpEvals1[i] <- 1/6*(delX1 + delX2)*(
							(2 - delX2/delX1)*y0 +
							(delX1 + delX2)^2/(delX1*delX2)*y1 +
							(2 - delX1/delX2)*y2)
	}

	
	#calculates simpson's rule beginning with the second node and 
	#leaving off the last one
	pieces2 <- pieces - 1
	nodes2 <- nodes[c(-1, -len)]
	values2 <- values[c(-1, -len)]
	
	simpEvals2 <- vector(mode = "numeric", length = pieces2)

	for(i in seq(1,pieces2))
	{
		idx0 <- 2*i - 1
		idx1 <- 2*i
		idx2 <- 2*i + 1
		
		y0 <- values2[idx0]
		y1 <- values2[idx1]
		y2 <- values2[idx2]
		
		delX1 <- nodes2[idx1] - nodes2[idx0]
		delX2 <- nodes2[idx2] - nodes2[idx1]	
	
		simpEvals2[i] <- 1/6*(delX1 + delX2)*(
							(2 - delX2/delX1)*y0 +
							(delX1 + delX2)^2/(delX1*delX2)*y1 +
							(2 - delX1/delX2)*y2)
	}
	
	
	#performs trapezoidal rule on first interval
	trap1 <- 1/2*(values[1] + values[2])*(nodes[2] - nodes[1])

	
	#gets running total
	cumInt <- vector(mode = "numeric", length = len)
	cumInt[1] <- 0
	cumInt[2] <- trap1
	cumInt[seq(3,len,2)] <- cumsum(simpEvals1)
	cumInt[seq(4,(len-1),2)] <- cumsum(simpEvals2) + trap1
	
	return(cumInt)
	
	
}
	
	