For example
ln(e)=(1-rhoe)*ln(steady_state(e))+rhoe*e(-1)+epsilone;
ln(z)=(1-rhoz)*ln(steady_state(z))+rhoz*z(-1)+epsilonz;
prevents computing the steady state numerically as you are trying to endogenously compute a steady state that needs to be specified by you. Also, you specified initial values of 0 for e,z, but you are taking their logs.
Also, see An infinity of steady states with Taylor rules