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

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