You are probably going to need to work at the MuPAD level, even if only by calling evalin(symengine) or feval(symengine)
However, you need to understand that expression trees need to match exactly to do this kind of manipulation. You can might be able to find all the subtrees that are 1/(ti-tj)+1/(ti-tk) and replace them, but 2/(ti-tj)+2/(ti-tk) is a different subtree that is not composed of the subtree _mult(2, 1/(ti-tj)+1/(ti-tk)) . Likewise, if you were to be searching for A*x^5 + B*x^4 + C*x^3 + 1 then the tree A*x^5 + B*x^4 + C*x^3 + 2 will not match it because although we as humans can easily see that is (A*x^5 + B*x^4 + C*x^3 + 1) + 1, MuPAD is not going to know that. And although you can search for a pattern A*x^5 + B*x^4 + C*x^3 + A_NUMERIC_CONSTANT, as humans we would also like to be able to automatically match A*x^5 + B*x^4 + 2*C*x^3 + 1 as being (A*x^5 + B*x^4 + C*x^3 + 1) + C*x^3 .. and then you want A*x^5 + B*x^4 + 2*C*x^3 + 5 to be recognized as (A*x^5 + B*x^4 + C*x^3 + 1) + (C*x^3 + 4)... when do you stop in the effort to match?
You might also find that in some of the subtrees, the same values are present but in a different order, depending on exactly how the trees got constructed. Subtrees are not sorted automatically (especially multivariate ones), so subtree x + y + 1 might be y + 1 + x . Negatives can get messy, -x-y might be _plus(_mult(-1,x),_mult(-1,y)) or might be _mult(-1,_plus(x,y)) or what you observe visually as (-x-y)/z might be stored as _mult(_plus(x,y),_mult(-1,_power(y,-1)))
