Working with a single equation, split up into several

General Usage
2

I kind-of solved the problem with ugly hacks. If anyone is interested, here is my haphazard solution, which just calls Symbolics.substitute until all the helper variables are gone.

using Symbolics

"""
Substitutes expressions of `target_eqn` as specified by equations `source_eqns`.
For each equation in `source_eqns`, substitutes its `lhs` by its `rhs` until none of the `lhs` are left in `target_eqn`.

TODO this is just a crappy prototype!
"""
function substitute_all(target_eqn, source_eqns,)
    all_eqns = [target_eqn, source_eqns...]

    eliminated_symbols = Set(map(x->x.lhs, source_eqns))
    length(eliminated_symbols) == length(source_eqns) || error("Redundant equations!")

    for i in 2:length(all_eqns)
        all_eqns = substitute.(all_eqns, (Dict(all_eqns[i].lhs => all_eqns[i].rhs),))
    end

    return all_eqns[1]
end


# test the function using above example

@variables u v w x y z a b

eqn_to_recover = z ~ a * u^2 + b * (v + 2w + 3x + 4y)

# define split-up version
helper_variables = @variables h_1, h_2, h_3

target_eqn = z ~ h_1 + b * (h_2 + h_3)

source_eqns = [
    h_1 ~ a * u^2
    h_2 ~ v + 2w
    h_3 ~ 3x + 4y
]

substitute_all(target_eqn, source_eqns,) == eqn_to_recover || error("Failure")