syms x y real           %criar var simbolicas
system=[-y - x^3;x^5]     %sistema
eqlb=solve(system);
eqlb=[eqlb.x eqlb.y]

V=x^6/3+y^2 

%dV = dV/dr*r_ponto + dV/dy*y_ponto
dV=gradient(V,[x y])'*system;
dV=expand(dV)

ezsurf(dV,[-0.5 0.5],[-0.5 0.5])