Tangent Line Between Two Functions

Question

Given two functions I want to find the tangent line(s) between them. The formula/system for this, if given two arbitrary functions f(x) and g(x), is:

f(x)-x*f'(x)-g(y)+y*g'(y)=0
f'(x)-g'(y)=0

where when solved, say x=a and y=b, gives two points: (a,f(a)) and (b,g(b)), such that a line can be drawn between them. The functions I want to find the tangent line between are x^-x and its first derivative (x^-x(-ln(x)-1)). The code I have is:

import sympy as sp
import numpy as np
import scipy.optimize as opt

n=1
s1,s2=sp.symbols('s1 s2')

def f1(z1):
return z1**-z1


def myFunction(z):
   x=z[0]
   y=z[1]
   q1=f1(s1)-s1*sp.diff(f1(s1),s1,n)- 
   sp.diff(f1(s2),s2,n)+s2*sp.diff(f1(s2),s2,n+1)
   q2=sp.diff(f1(s1),s1,n)-sp.diff(f1(s2),s2,n+1)

   F=np.empty((2))
   F[0]=abs(q1.subs([(s1,x),(s2,y)]))
   F[1]=abs(q2.subs([(s1,x),(s2,y)]))
   return F
zGuess=np.array([0.6,1.3])
z=opt.fsolve(myFunction, zGuess)
print(z)

The solution should be x=0.59515 and y=1.37866 but it says "Can't calculate 1st derivative wrt 0.6". I don't know how to have it input the initial values for solving the system after taking the derivative with respect to x and y respectively. How would I get it to solve for these values while being able to change n so that the tangent line can also be found between the first derivative and second derivative and so on? Here is a picture showing what it looks like.

Answer 1

I figured out how to solve this problem. Here is the solution below.

import sympy as sp
import numpy as np
import scipy.optimize as opt
s1,s2=sp.symbols('s1 s2')

def f1(z1):
    return z1**-z1

def myFunction(z,n):
    x=z[0]
    y=z[1]

    if n % 2 ==0:
        w1=1
        w2=-1

    else:
        w1=-1
        w2=1
    q1=w2*sp.diff(f1(s1),s1,n-1)-w2*s1*sp.diff(f1(s1),s1,n)- 
w1*sp.diff(f1(s2),s2,n)+w1*s2*sp.diff(f1(s2),s2,n+1)
    q2=w2*sp.diff(f1(s1),s1,n)-w1*sp.diff(f1(s2),s2,n+1)

    F=np.empty((2))
    F[0]=abs(q1.subs([(s1,x),(s2,y)]))
    F[1]=abs(q2.subs([(s1,x),(s2,y)]))
    return F


j1=1
j2=10
z=np.empty((j2+1,2))
for i in range(j1,j2+1):
    n=i
    a=float(0.5211*n+0.184)
    b=float(0.916745*a+0.794333)
    zGuess=np.array([a,b])
    z[i,:]=opt.fsolve(myFunction, zGuess, n, maxfev = 6000)
    print(z[i,:])

This solution gives the first ten pairs of points. The output when done is

[ 0.59515445  1.37866112]
[ 1.19255583  1.87697191]
[ 1.75662656  2.37914507]
[ 2.29484826  2.87469058]
[ 2.81278951  3.36119549]
[ 3.31426591  3.83840029]
[ 3.80198665  4.30671276]
[ 4.277955    4.76673553]
[ 4.74370426  5.21909897]
[ 5.20044161  5.66440097]