ideal membership and solution

edit

Yes, this can be done algorithmically. First you compute a Groebner basis $\langle g_1, \ldots, g_r \rangle$ of $I$, and remember the expressions $g_i = a_i f + b_i g$. Then you do multivariate polynomial division of $h$ by $\langle g_1, \ldots, g_r \rangle$. This yields $$h = c_1 g_1 + \ldots + c_r g_r = (c_1a_1 + \ldots c_ra_r)f + (c_1b_1 + \ldots c_rb_r)g.$$

In your example one Groebner basis of $I$ is $\langle f,g,f_1,f_2,f_3 \rangle$ where $f_1 = S(f,g) = f-zg$, $f_2 = S(f_1,f) - zg = (y-z)f - z(y+1)g$, $f_3 = S(f_1,f_2) = (x - z(y-z))f + z(z(y+1) - x)g$. Multivariate polynomial division of $h$ by $\langle f,g,f_1,f_2,f_3 \rangle$ yields $$h = f + f_2 = (y-z+1)f - z(y+1)g.$$

This can all easily be automated in Sage (but I don't have time to do it now).