# \hat u(k) computed using Newton algorithm
function hatu_newton(order,rho,a0,d,v,maxiter,tolerance)

  # compute gaussian quadrature weights
  weights=gausslegendre(order)

  # initialize V and Eta
  (V,V0,Eta,Eta0)=init_veta(weights,d,v)

  # initialize u, V and Eta
  u=zeros(Complex{Float64},order)
  for j in 1:order
    # transformed k
    k=(1-weights[1][j])/(1+weights[1][j])
    if d==2
      u[j]=2*pi*a0*rho/k
    elseif d==3
      u[j]=4*pi*a0*rho/k^2
    end
  end

  # iterate
  for i in 1:maxiter-1
    new=u-inv(dXi(u,V,V0,Eta,Eta0,weights,rho,d))*Xi(u,V,V0,Eta,Eta0,weights,rho,d)

    if(norm(new-u)/norm(u)<tolerance)
      u=new
      break
    end
    u=new
  end

  return (u,en(u,V0,Eta0,rho,weights,d)*rho/2)
end

# Xi(u)=0 is equivalent to the linear equation
function Xi(u,V,V0,Eta,Eta0,weights,rho,d)
  order=length(weights[1])

  # init
  out=zeros(Complex{Float64},order)

  # compute E before running the loop
  E=en(u,V0,Eta0,rho,weights,d)

  for i in 1:order
    # k_i
    k=(1-weights[1][i])/(1+weights[1][i])
    # S_i
    S=V[i]-1/(rho*(2*pi)^d)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights)
    # X_i
    X=k^2/(2*E*rho)

    # U_i
    out[i]=u[i]-S/(2*E*(X+1))*Phi(S/(E*(X+1)^2))
  end

  return out
end

# derivative of Xi (for Newton)
function dXi(u,V,V0,Eta,Eta0,weights,rho,d)
  order=length(weights[1])

  # init
  out=zeros(Complex{Float64},order,order)

  # compute E before the loop
  E=en(u,V0,Eta0,rho,weights,d)

  for i in 1:order
    # k_i
    k=(1-weights[1][i])/(1+weights[1][i])
    # S_i
    S=V[i]-1/(rho*(2*pi)^d)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights)
    # X_i
    X=k^2/(2*E*rho)

    for j in 1:order
      y=(weights[1][j]+1)/2
      dS=-1/rho*(1-y)*Eta[i][j]/(2*(2*pi)^d*y^3)*weights[2][j]
      dE=-1/rho*(1-y)*Eta0[j]/(2*(2*pi)^d*y^3)*weights[2][j]
      dX=-k^2/(2*E^2*rho)*dE
      out[i,j]=(i==j ? 1 : 0)-(dS-S*dE/E-S*dX/(X+1))/(2*E*(X+1))*Phi(S/(E*(X+1)^2))-S/(2*E^2*(X+1)^3)*(dS-S*dE/E-2*S*dX/(X+1))*dPhi(S/(E*(X+1)^2))
    end
  end

  return out
end

# energy
function en(u,V0,Eta0,rho,weights,d)
  return V0-1/(rho*(2*pi)^d)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*u,0,1,weights)
end