SymPy `subs` Not Doing Anything

I have a differential equation for which I use sympy.solvers.ode.dsolve to solve, I get out

                ___________           ___________
           -x⋅╲╱ E - V_max        x⋅╲╱ E - V_max 
ψ(x) = C₁⋅ℯ                 + C₂⋅ℯ               

From (I put some of the code to I used to generate this equation at the end):

                         2          
                        d           
-ψ(x)⋅E + ψ(x)⋅V_max + ───(ψ(x)) = 0
                         2          
                       dx           

This is all well and good, the problem comes when I know that C₁ and C₂ happen to be equal and want to substitute one for the other. So I try something like

psi_high.subs( sp.Symbol( "C_2" ), sp.Symbol( "C_1" ) )

However it just comes out the same as before

                ___________           ___________
           -x⋅╲╱ E - V_max        x⋅╲╱ E - V_max 
ψ(x) = C₁⋅ℯ                 + C₂⋅ℯ               

I am thinking this may be a memory issue, that a reference to a sympy.Symbol object must refer to not only a sympy.Symbol object with the same value/symbol but which also must be the same underlying object.

This is only speculation (but I can say, do psi_high.subs( x, 0 ) and it works), but my question is how do I resolve it?

Curiously, it seems to work here (I did try this using the sympy.symbols function and by enclosing the symbol references in a tuple and list like shown in the question)

Thanks!

well_length = sq.Quantity( 'L' )
highest_potential = sq.Quantity( "V_max" )
x = sp.Symbol( 'x' )
m = sq.Quantity( 'm' )
hbar = sq.Quantity( "hbar" )
total_energy = sq.Quantity( 'E' )
inverse_total_energy = 1.0 / total_energy
psi_symbol = ud.lookup( "GREEK SMALL LETTER PSI" )
psi = sp.Function( "psi" )
second_derivative = sp.Derivative( psi( x ), x, 2 )
make_shrodinger_left = lambda potential, psi_parameter : ( second_derivative + ( psi( psi_parameter ) * potential ) )
make_shrodinger_right = lambda psi_parameter : total_energy * psi( psi_parameter )
make_psi_equal = lambda input_value, value : sp.Eq( psi( sp.Eq( x, input_value ) ), value )
set_equal = lambda to_set, value : sp.Eq( to_set, value )
shrodinger_left_high = sp.simplify( make_shrodinger_left( highest_potential, x ) )
shrodinger_right = make_shrodinger_right( x )
high_diff = sp.simplify( set_equal( shrodinger_left_high - shrodinger_right, 0 ) )

>Solution :

Here’s a simpler example:

In [3]: eq = Eq(f(x).diff(x, 2), 0)

In [4]: eq
Out[4]: 
  2          
 d           
───(f(x)) = 0
  2          
dx           

In [5]: sol = dsolve(eq)

In [7]: sol
Out[7]: f(x) = C₁ + C₂⋅x

We can inspect these symbols:

In [8]: sol.free_symbols
Out[8]: {C₁, C₂, x}

In [9]: [s.name for s in sol.free_symbols]
Out[9]: ['C2', 'C1', 'x']

Note that there are no underscores in the symbol names. What we want to do then is:

In [10]: sol.subs(Symbol("C1"), Symbol("C2"))
Out[10]: f(x) = C₂⋅x + C₂