Multi qubit expectation values

Please visit the following post for more clarification:

Explanation post and code

A dummy code:

def my_quantum_function(x, y):
    qml.RZ(x, wires=0)
    qml.CNOT(wires=[0, 1])
    qml.RY(y, wires=1)
    return [qml.expval(qml.PauliZ(i)) for i in range(2)]

Here, how are the two expvals are calculated theoretically?