How do I perform a change of variables for a pde

edit

With some insights from the following posts:

substitute expression instead of formal function symbol

How do I get an ordered list of a symbolic functions arguments?

Call a function with argument list in python (non-sage link)

I've come up with a function, transform_func_x(expr,f,rels), to change the dependent variables of an expression. Below is the function itself and below that how I've used it (all copy and pasted from a sage notebook). Any suggested improvements are welcome.

def transform_func_x(expr,f,rels):
    r"""
    Transform the dependent variables of a symbolic function

    INPUT:

    - ``expr`` - symbolic function.  Usually a relational expression
      with differential operators such as a partial differential equation.

    - ``f`` -  expression defining a function (e.g. previously defined
      as f = function ('u', z, t))

    - ``rels`` - list of relational expressions giving variable
      transformations.  For example to change y(z, t) to y(Z, T_v), rels
      could be [Z == z/H, T_v = c_v * t / H^2].  To change y(z, t) to y(Z, t),
      rels could be [Z == z / H]

    .. note::

       1. Transformations will not be performed for any relationship where the
       LHS is already one of the variables in f. For, example if
       f = u(Z, t) and rels = [Z == z / H, T_v == c_v * t / H^2] then the
       "Z == z / H" transformation will be ignored.
       2. After applying this function to expr, you will probably want to 
       update f by applying transform_func_x to f itself 
       i.e transform_func_x(f,f,rels)
    """
#Limit#all#lines#to#a#maximum#of#79#characters#################################
    y = f.operator()
    x = f.operands()
    nx = f.number_of_operands()

    new_arg_list = x[:]
    """used to transform variables within f to functions
       of those variables. e.g. u(z, t) --> u(z / H, t)"""
    subs_dict1 = {}
    """used to change name of transformed variables within f.
       e.g. u(z / H, t) --> u(Z, t)"""
    subs_dict2 = {}
    """used to transform variables outside of f.
       e.g. z *u(Z, t)--> z * H * u(Z, t)"""
    for v in rels:
        if v.left() not in x:
            """i.e. exclude transformations where LHS of
               rels is already a variable of f"""
            for i in range(nx):
                # find which variable the transformation applies to
                if x[i] in v.right().variables():
                    xi=i
                    break
            new_arg_list[xi] = v.right()
            subs_dict1[v.right()]=v.left()
            subs_dict2[x[xi]]=solve(v,x[xi],solution_dict=true)[_sage_const_0 ][x[xi]]

    # make dummy function for transforming variables of f
    dummy_arg_list = [var("x{0:d}".format(i)) for i in range(nx)]
    dummy_y(*dummy_arg_list) = function(str(y),*new_arg_list))

    return (expr.substitute_function(y,dummy_y).subs(subs_dict1).
            subs(subs_dict2))

reset
var('z,t,H,T_v,Z,c_v,k_v,m_v, g_w')
f=function('u',z,t)
pde = k_v*f.diff(z,z)/g_w == m_v*f.diff(t)

print('Variable transformations')

print('PDE with u(z,t):')
rels=[]
eq1=transform_func_x(pde,f,rels)
f=transform_func_x(f,f,rels)
print(eq1)
print('Term to isolate:')
myterm = f.diff(f.operands()[0],2)
print (myterm)
print('Coefficient of term on LHS:')
mycoeff=eq1.left_hand_side().coeff(myterm)
print(mycoeff ...