function [f,dfdx,dfdy] = cubic(x,y)

f = y - (x+1).*x.*(x-1);
dfdx = -3.*(x.^2) + 1;
dfdy = 1;