This is the python script used to do symbolic calculations of derivatives for article 79649

#!/usr/bin/env python3

# ================================================================
#  This script does symbolic calculations for article 79649
#
#  Written by Randy S (2022)
# ================================================================

#import numpy as np
import copy

which_equation = 2
# which_equation options:
#  1 = the n-th derivative of V with respect to \phi
#  2 = \beta_n

# ===================================================================
def show_term_1( term ):

    output = '&'

    if term['coef'] > 0:
        output += '+'
    output += str( term['coef'] ) + '\\frac{'
    Vexp = 1
    for k in range(len(term['Vfactors'])):
        if Vexp > 1:
            Vexp -= 1
        else:
            output += 'V_{' + str( term['Vfactors'][k] ) + '}'
            Vexp = 1
            while k+Vexp < len(term['Vfactors']) and term['Vfactors'][k+Vexp] == term['Vfactors'][k]:
                Vexp += 1
            if Vexp > 1:
                output += '^{' + str(Vexp) + '}'
    output += '}{D^{' + str( term['Dexp'] ) + '}}'
    print( output )
    print( '\cr')

# ===================================================================
def show_terms_1( n, terms ):

    print( '\\begin{align*}' )
    print( '\pl_\phi^{' + str(n) + '}\log(D) &= ' )
    print( '\cr')
    for k in range(len(terms)):
        show_term_1( terms[k] )
    print( '\end{align*}' )

# ===================================================================
def show_term_2( term ):

    any_odd = False

    output = '&'

    if term['coef'] > 0:
        output += '+'
    output += str( term['coef'] ) + '\\frac{'
    Vexp = 1
    for k in range(len(term['Vfactors'])):
        if Vexp > 1:
            Vexp -= 1
        else:
            if term['Vfactors'][k] % 2 == 1:
                any_odd = True
            output += 'g_{' + str( term['Vfactors'][k] ) + '}'
            Vexp = 1
            while k+Vexp < len(term['Vfactors']) and term['Vfactors'][k+Vexp] == term['Vfactors'][k]:
                Vexp += 1
            if Vexp > 1:
                output += '^{' + str(Vexp) + '}'
    output += '}{(1+g_2)^{' + str( term['Dexp'] ) + '}}'

    if not any_odd: # only show terms without any factors of g_k with odd k
        print( output )
        print( '\cr')

# ===================================================================
def show_terms_2( n, terms ):

    if n%2 == 1:
        return  # only show results for n even

    print( '\\begin{align*}' )
    print( '\\beta_{' + str(n) + '} &= ' )
    print( '(' + str(n) + ' - ' + str(int(n/2-1)) + 'd)g_{' + str(n) + '} ' )
    print( '\cr')
    for k in range(len(terms)):
        show_term_2( terms[k] )
    print( '\end{align*}' )

# ===================================================================
def show_terms( n, terms ):

    if which_equation == 1:
        show_terms_1( n, terms )
    else:
        show_terms_2( n, terms )

# ===================================================================
def take_deriv_of_one_term( term ):

    result = []

    # append the term that comes from the derivative hitting the denominator
    newfactors = copy.deepcopy( term['Vfactors'] )
    newfactors.append(3)
    result.append( {'Dexp': term['Dexp']+1, \
                    'Vfactors': newfactors, \
                    'coef': -term['Dexp'] * term['coef'] } )

    # append the terms that come from the derivative hitting each factor in the numerator
    for k in range(len(term['Vfactors'])):
        newterm = copy.deepcopy( term )
        newterm['Vfactors'][k] = newterm['Vfactors'][k] + 1
        result.append( newterm )

    return result

# ===================================================================
def consolidate_result( result, result_single ):

    for k in range(len(result_single)):
        result_single[k]['Vfactors'].sort( reverse = True )

    for k in range(len(result_single)):
        found_match = False
        for j in range(len(result)):

            if result_single[k]['Vfactors'] == result[j]['Vfactors']:
                # these terms have the same factors, so consolidate them
                result[j]['coef'] += result_single[k]['coef']
                found_match = True
        if not found_match:
            result.append( result_single[k] )

# ===================================================================
def take_deriv( terms ):

    result = []
    for k in range(len(terms)):
        result_single = take_deriv_of_one_term( terms[k] )
        consolidate_result( result, result_single )

    return result

# ===================================================================
#   Main program

# start with the first derivative of log(D)
terms = [{'Dexp': 1, 'coef': 1, 'Vfactors': [3]}]
show_terms( 1, terms )

# calculate higher derivatives iteratively
for n in range(2,11):
    terms = copy.deepcopy( take_deriv( terms ) )
    show_terms( n, terms )