How to find a pair of divisors as close as possible to each other?

Note that although there could be more efficient algorithms, what you found is not a bad algorithm. I will analyze this a little bit.

Let us write down the problem formally as a constrained discrete optimization problem: $$ (x^star, y^star) = argmin_{(x,y) in mathbb{Z_{+}}} Big(E triangleq big(y-xbig)^2 + lambda(x*y - n)^2 Big) $$ Note that for positive integers the $L_2$ norm will attain the same minimum as $L_1$. So I'll stick to it. $lambda$ is a Lagrange multiplier and the right-most term is the non-linear constraint. Even if we can compute the partial derivatives and use the method of Lagrange multipliers, the problem is not easy to solve due to the parameters restricted to a discrete space. Hence we will relax the integer constraint to real numbers and solve the problem for $mathbb{R}$. The relaxed optimum would be attained as follows: $$ x^star = argmin_{x in mathbb{R}} Big(tilde{E} triangleq big(x-frac{n}{x}big)^2 Big) $$ We achieved this by substitution the constraint : $y gets frac{n}{x}$. Note that this is also suggested by the Largange optimization: $$ frac{partial E}{partial lambda} = xy - n^2 = 0 implies x=frac{n}{y} text{ or } y=frac{n}{x}. $$

Now differentiating w.r.t. $x$: $$ begin{align} frac{d tilde{E}}{dx} = 2big(x-frac{n}{x}big) &= 0 x^2 &= n x = sqrt{n} end{align} $$ gives the optimal solution for the relaxed problem. For negative numbers $x = -sqrt{n}$ is also the optimum. This we could already detect from the sign of the input. Having solved the relaxed problem, a typical optimizer would operate by projecting the solution into the feasible set and/or applying a projected gradient descent. The initial feasible solution is $x_0 = lfloor sqrt{n} rceil$. $lfloor cdot rceil$ is the rounding operation or in fancy words the projection onto the integers. From here, it is possible to walk on the Euclidean discrete manifold of 1-dimensional integers by your method: taking steps of $+1$ or $-1$ depending on the direction. That would amount to a gradient descent by the step size $h=1$.

And these steps are what you are following implicitly when arriving at the algorithm you proposed. Unfortunately, the worst case complexity of this algorithm is attained when the number is prime and we have only 2 divisors: $1$ and $n$, the number itself. Therefore, the worst case complexity reaches to $O(sqrt{n})$. As $nrightarrow infty$, this would resemble a linear search (it is worse than $log$). Now, the set of integers is a totally ordered set. Using this property, instead of using a linear search, one could think of a binary search variant. But as clarified by @WolfgangBangerth, such an algorithm would not exist bounding the complexity in: $O(log(sqrt{n}))<Cleq O(sqrt{n})$.

Even though this might not be the exact answer you were looking for, I thought it would give a perspective to develop better solutions. For instance, the integer programming community might have a better solution.