Matlab - sistema ecuaciones

Hola a todos, me gustaría resolver un sistema de 3 ecuaciones con 3 incógnitas. El sistema podría ser el siguiente:

x= y+A*y^2*z ec. 1
y=z*x^2-B*x ec.2
z*x-z^2*y+C*x^2*y^2*x=0-------ecuación de la que no puedo despejar z ec.3

Mi idea es hacer lo siguiente:
Supongo z (Ahora tengo un sistema de 2 ecuaciones: ec1, ec2 con 2 incógnitas)
Supongo x
Con la ecuación 2 calculo y
Con la ecuación 1 calculo x
Compruebo x
Compruebo z

Todo esto lo haría con bucles, pero me parece que estoy complicando mucho el problema y que tendría una solución más fácil, ¿Me podeis ayudar?

Muchas gracias

Para complicar más las cosas, el sistema ya no es de 3 ecuaciones sino de 4. Las incógnitas son k2,H,A,I.

ecuación 1: función de k2,H,A,I
ecuación 2: función de k2,H,A,I
ecuación 3: H = función de A,I
ecuación 4 I= función de A,H

no puedo simplificar ninguna ecuación más, ¿como lo resolveríais?

El sistema de ecuaciones que os había comentado es el que está en esta página:

http://rinconmatematico.com/foros/index.php?topic=20626.new#new

Hola, estoy intentando resolver el sistema, mi idea es la siguiente:

1. Defino lambda y epsilon porque son constantes

2. Defino a (cómo está contenida entre cero y uno, le doy un valor muy pequeño, 0.001)

3. Defino h (cómo es mayor que uno, le doy de valor 1.001)

4. Escribo la ecuación (3), es decir, la función que define a "h"

5. Defino k (cómo puede ser negativa, le doy de valor -2)

6. Escribo la ecuación (2) en la que se define a "k"

7. Le digo que si la ecuación(3) es distinta que "h" que coga cómo valor de h = h+0.001

8. Lo mismo para k

Os escribo lo que he hecho:

lambda=1;
epsilon=0.001;
a=0.001;
h=1.001;
h1=1+(1-a)*(1+a-(1+h)*lambda*tand(epsilon))/(lambda((1+h)*lambda+(1+a)*tand(epsilon)));
k=-2;
k1=(1+tand(epsilon)^2)*(-lambda^2*(1-a^2)-k(1-a^2)/(lambda(1+h)^2))/(1+(1+a)/(lambda(1+h))*tand(epsilon))^2;
while abs(h1-h)>0.001
h=h+0.001
end
while abs(k1-k)>0.001
k=k+0.001
end
display(a,h1,k1)

Me sale de error:
Subscript indices must either be real positive integers or logicals.

Error in ==> sistema at 5
h1=1+(1-a)*(1+a-(1+h)*lambda*tand(epsilon))/(lambda((1+h)*lambda+(1+a)*tand(epsilon)));

¿Alguien sabe que significa?

Muchas gracias