function xx_sol=Resuelve_con_Newton(niter,EPSILON,eps,...
    xx,ns,beta,alpha,a,LT,KT,theta)

dff=zeros(size(1,2*ns+1));

% inicio del newton
for  iter=1:niter
    [ff]=ff_eq(xx,ns,beta,alpha,a,LT,KT,theta);
    diff=0;
    for i=1:ns*2+1
        diff=diff+abs(ff(i));
    end

    if diff <= EPSILON
        disp('Enhorabuena hemos convergido');
        break
    end
    % computamos el proximo punto
    %eps=10^(-5);
    for j=1:2*ns+1
        xx(j)=xx(j)+eps;
        [ff1]=ff_eq(xx,ns,beta,alpha,a,LT,KT,theta);
        for i=1:2*ns+1
            dff(i,j)=(ff1(i)-ff(i))/eps;
        end
        xx(j)=xx(j)-eps;
    end

    if iter <=5
        xscal=.5;
    else
        xscal=1;
    end

    dffinv=dff^(-1);
    for j=1:2*ns+1
        for i=1:2*ns+1
            xx(i)=xx(i)-xscal*dffinv(i,j)*ff(j);
        end
    end
end

xx_sol=xx;
