This algorithm finds the solution of an non-linear equation. The equation is the catenary, which must pass the point (x0, y0). There are 0, 1 or 2 solutions for this problem

findcatenary.py

######################################
#             Description            #
######################################
#
#  This algorithm finds the solution of an non-linear equation
#  The equation is the catenary, that is like
#            h(x) = a * cosh(x/a)
#  We want to find the value of "a" such that h(x0) = y0
#  Where (x0, y0) is known
#
#  In this problem, there's four types of solution
#     if x0 == 0:  one solution
#     elif y0/x0 < ftstar:  zero solutions
#     elif y0/x0 == ftstar:  one solution  # almost impossible to get this case
#     elif y0/x0 > ftstar:  two solutions
#  The value of ftstar is constant, approx 1.50888. See the link math.stackexchange.com below 
#
#  Link to math.stackexchange
#  https://math.stackexchange.com/questions/4434450/find-parameter-to-catenary-interpolate-a-specific-point/4434815#4434815

from scipy.special import lambertw
import numpy as np
from matplotlib import pyplot as plt
from numpy import linalg as la

sinh = np.sinh
cosh = np.cosh

def f(t):
    return np.cosh(t) / t

def df(t):
    return np.sinh(t) / t - f(t) / t

def g(t, k):
    """
    This function is to iterate, using Newton's method to find the solution
    t_{n+1} = t_{n} - f(t_{n})/f'(t_{n})
    t_{n+1} = t_{n} - g(t_{n})
    ----> g(t) = f(t)/f'(t)
    """
    return t * (np.cosh(t) - k * t) / (t * np.sinh(t) - np.cosh(t))


x0, y0 = 1, 2
tstar = 1.1996786402577338339
fstar = f(tstar)
print("tstar = ", tstar)
print("fstar = ", fstar)

if x0 == 0:
    a1 = y0
    a2 = None
    print("a = ", a1)
else:
    k = y0 / x0
    print("    k = ", k)
    if k < fstar:
        msg = "k is too low. Impossible find catenary to solve it\n"
        msg += "k must be > %.4f" % fstar
        print(msg)
        a1 = None
        a2 = None
    else:
        t1 = 1 / k + 1 / (2 * k**3) + 13 / (24 * k**5) + 541 / \
            (720 * k**7) + 9509 / (8064 * k**9)
        print("frist estimate t1 = ", t1)  # Frist estimate. Now we use Newton
        for w in range(5):  # Only 5 iterations
            t1 -= g(t1, k)
        print(" last estimate t1 = ", t1)
        a1 = x0 / t1
        print("a1 = ", a1)
        print("")

        t2 = -lambertw(-1 / (2 * k), -1)
        if la.norm(t2 - np.real(t2)) < 1e-8:
            t2 = np.real(t2)  # Frist estimate. Now we use Newton
            print("frist estimate t2 = ", t2)
            for w in range(5):  # Only 5 iterations
                t2 -= g(t2, k)
            print(" last estimate t2 = ", t2)
            a2 = x0 / t2
            print("a2 = ", a2)
        else:
            a2 = None

##############################################
#                   PLOTTING                 #
##############################################

ymin = 0
ymax = np.abs(y0)
if a1 is not None:
    ymax = np.max([ymax, np.abs(a1)])
if a2 is not None:
    ymax = np.max([ymax, np.abs(a2)])
ymax *= 1.5
L = np.log(2 * ymax)
L = np.max([L, 1.5 * x0])
xplot = np.linspace(-L, L, 129)
plt.scatter(x0, y0)
if a1 is not None:
    plt.plot(xplot, a1 * cosh(xplot / a1), label="Frist solution")
if a2 is not None:
    plt.plot(xplot, a2 * cosh(xplot / a2), label="Second solution")
plt.plot((-L, 0, L), (fstar * L, 0, fstar * L),
         ls="dashed", label="Limit line")
plt.ylim([0, ymax])
plt.xlim([-L, L])
plt.legend()
plt.show()