Is a function in an ideal? Verification by hand and Macaulay 2 Suppose 
$$f_1=-4x^4y^2z^2+y^6+3z^5,$$ 
$$f_2=-4x^2y^2z^2+y^6+3z^5,$$
$$f_3=4x^4y^2z^2+y^6+3z^5,$$ 
$$f_4=4x^2y^2z^2+y^6+3z^5$$
and
$$I=\langle xz-y^2,x^3-z^2\rangle\subset\mathbb C[x,y,z].$$ 
Is $f_i\in I?$

The answer is Yes in some cases.
The question can be checked with Macaulay 2: when the remainder is zero with respect to the Gröbner basis like (R=QQ[x,y,z]; fi=...; I=ideal(x*z-y^2,x^3-z^2); G=gb(I);f%G returning zero, $f_i\not\in I$. 
Division with respect to the elements in ideal or with respect to polynomials in the basis does not result into the case with zero remainder as demonstrated in the trials with Lex, RevLex and with polynomials in Gröbner basis and with respect to polynomials in the ideal (in larger resolution here, here and here) where all trials fail because not finding the zero remainder case.

Checking whether a polynomial is in an ideal can be done by hand such as manually with division algorithm and with Macaulay2 such as finding the correct factors, decomposition of the polynomial with quotients and GR elements. So
How can you check that a polynomial is in an ideal by hand and computationally such as Macaulay 2?
P.s. Other questions on learning basic computational algebraic geometry in chat.
 A: I tried the lex order. In your computation you misinput $-4x^2y^2z^2+y^6+3z^5$ as $4x^2y^2z^2+y^6+3z^5$. Otherwise using lex order can easily give you remainder $0$. 
$$
\require{enclose}
\begin{array}{cc}
                    &-4z^2 + 1  \\[-3pt]
x^3-z^2& \enclose{longdiv}{-4x^2y^2z^2+y^6+3z^5} \\[-3pt]
x^2y^2-z^3 &        \underline{-4x^2y^2z^2+4z^5}\phantom{+3z^5} \\[-3pt]
xy^4-z^4 &                    y^6-z^5  \\[-3pt]
xz-y^2 &         \underline{y^6-z^5}\\[-3pt]
y^6-z^5& 0
\end{array}
$$
Apparently then $f=-4z^2(x^2y^2-z^3)+(y^6-z^5)$.
Here is a program in Sage that I adopted from: http://www.math.utah.edu/~carlson/cimat/lecture1.pdf:
def div(f,g):
n = len(g)
p, r, a = f, 0,  [0 for x in range(0,n)]
while p != 0:
    i, divisionoccured = 0, False
    while i < n and divisionoccured == False:
        if g[i].lt().divides(p.lt()):
            a[i] = a[i] + p.lt()//g[i].lt()
            p = p - (p.lt()//g[i].lt())*g[i]
            divisionoccured = True
        else:
            i = i + 1
    if divisionoccured == False:
        r = r + p.lt()
        p = p - p.lt()
return a, r

Results from running the above algorithm:
P.<x,y,z> = PolynomialRing(QQ,order='lex')
I = ideal(x*z-y^2,x^3-z^2)
B = I.groebner_basis(); B
f = -4*x^2*y^2*z^2+y^6+3*z^5
div(f,B)
[x^3 - z^2, x^2*y^2 - z^3, x*y^4 - z^4, x*z - y^2, y^6 - z^5]
([0, -4*z^2, 0, 0, 1], 0)

A: It is interesting that Macaulay2 solution in terms of ideal quotients differs from the solution by hand. Consider $f_2=-4x^2y^2z^2+y^6+3z^5$. M2 in the polynomial ring with Lex order provided 
$$f_2^{lex,M2}=(3x^2z^2-xy^2z-y^4)(xz-y^2)+(-3z^3)(x^3-z^2)$$
while by hand in inverse lexicographic order
$$f_2^{invlex}=(-4xy^2z-4y^4)(xz-y^2)-3(y^6-z^5)$$
and by KittyL in Lex order
$$f_2^{lex,Hand}=-4z^2(x^2y^2-z^3)+(y^6-z^5)$$
where $f_2^{lex,M2}$ should be equal to $f_2^{lex,Hand}$ but it is not: so M2 is using some different order in division and quotients. It is a task to investigate why M2 provides so different result.
Macaulay2
R=QQ[x,y,z,MonomialOrder=>Lex]; f=-4*x^2*y^2*z^2+y^6+3*z^5; 
I=ideal(x*z-y^2,x^3-z^2); G=gb(I);f//gens G
gens G
(3*x^2*z^2-x*y^2*z-y^4)*(x*z-y^2)+(-3*z^3)*(x^3-z^2)
f//gens G
matrix(gens G)*matrix(f//gens G)==(3*x^2*z^2-x*y^2*z-y^4)*(x*z-y^2)+(-3*z^3)*(x^3-z^2)


and computational aspect, more here, about finding the factors/coefficients for the polynomials so that the linear combination equals zero.
By hand
On KittyL's observations with corrections such as $f$, division algorithm with respect to inverse Lex on $f$ and dividing with the polynomials in the Gröbner basis of the ideal:

