Implementing Mathematica's solve and eliminate functions using SymPy

sympytools.py

import sympy as sp

def eliminate(system, symbols):
    """Eliminates given variables from a system of equations.
    
    Args:
        system: List of SymPy equations.
        symbols: List of SymPy symbols to eliminate. Must be shorter than list of equations.
            
    Returns:
        A list of SymPy equations with N_symbols fewer equations than the input. The new
        equations have had the given symbols eliminated.
    """
    new_system = system.copy()
    for symbol in symbols:
        solvedfor = None
        elimidx = None
        for idx, eq in enumerate(new_system):
            if symbol in eq.free_symbols:
                if solvedfor is None:
                    solvedfor = sp.solve(eq, symbol, dict=True)[0]
                    elimidx = idx
                else:
                    new_system[idx] = eq.subs(solvedfor)
        del new_system[elimidx]
    return new_system


def condense(system, symbols, elim):
    """Solves a system for a set of variables while eliminating others.
    
    Args:
        system: List of SymPy equations. Must have len(system) == len(symbols) + len(elim).
        symbols: List of SymPy symbols to solve for.
        elim: List of SymPy symbols to eliminate.
            
    Returns:
        A list of len(symbols) SymPy equations with the given variables eliminated.
        Equations have symbols to be solved for isolated on the left hand side.
    """
    return [sp.Eq(k, v) for (k, v) in sp.solve(eliminate(system, elim), symbols).items()]

This is nice. Are you still using this?