function [lambda,mu]=newton(G,JG,lambda0,mu0,l,beta,tol);
  X0=[lambda0;mu0];
  napaka=tol+1;
  while napaka>=tol
    X=X0-vrni_jg(X0,l,beta)\vrni_g(X0,l,beta);
    napaka=norm(X-X0);
    X0=X
  end
  lambda=X(1);mu=X(2);

