om1=10; phi1=pi/3;
om2=13; phi2=-pi/5; x0=[0;0];

k=1; x=x0; err=1e+10; tol=1e-10; kmax=10;
while ((k<=kmax) && (err > tol))
 k
 f  = [x(2)-sin(om1*x(1)+phi1);
     x(1)-sin(om2*x(2)+phi2)];
 J = [ -cos(om1*x(1)+phi1)*om1, 1;
      1, -cos(om2*x(2)+phi2)*om2 ];
 norm(f)
 dx = -J\f;
 err = norm(dx)
 x=x+dx;
 k=k+1;
endwhile