module roots
	interface
		subroutine findroot(x, dx, eps)
			real, dimension(:), intent(inout) :: x
			real, dimension(size(x)), intent(in) :: dx
			real, intent(in) :: eps
		end subroutine findroot
	end interface
end module


subroutine findroot(x, dx, eps)
	use matrix, only : qrdec, lineqnsystem, norm
	implicit none
	real, dimension(:), intent(inout) :: x
	real, dimension(size(x)), intent(in) :: dx
	real, intent(in) :: eps
	INTERFACE
		FUNCTION f(x)
			real, intent(in) ::x(:)
			real :: f(size(x))
		END FUNCTION
	END INTERFACE
	real :: lambda
	real, dimension(size(x)) :: z, Deltax, df, fx
	real, dimension(size(x),size(x)) :: J, Q, R
	integer :: n, i, k
	n = size(x)

	do
		fx = f(x)
	
		! Make Jacobi-matrix:
		J = 0.0
		do k = 1,n
			x(k) = x(k)+dx(k)
			df = f(x)-fx
			x(k) = x(k)-dx(k)
			do i = 1,n
				J(i,k) = df(i)/dx(k)
			enddo
		enddo
	
		! Find Deltax by solving the system J*Deltax=-f(x)
		call qrdec(J,Q,R)
		Deltax = lineqnsystem(Q,R,-fx)
	
		lambda = 1.0
		do
			z = x + lambda*Deltax
			lambda = lambda/2
			if (norm(f(z)) <= (1-lambda)*norm(fx) .or. lambda <= 0.01) exit
		enddo
		x = z
		
		if (norm(f(x)) < eps .or. norm(lambda*Deltax) < norm(dx)) exit
	enddo
end subroutine findroot