Symbolically evaluating jump conditions of expressions across interfaces

I am dealing with a problem in electromagnetism involving problems in different regions, and there is something that I haven’t been able to figure out yet. What I would like to be able to do is take an algebraic expression in terms of symbolic functions and constants, and give them labels depending on which side of an interface they are at (thus creating new functions/constants).

As a minimal working example, consider the quantity

$$frac{1}{mu}dfrac{partial B}{partial x} $$

where $mu$ is a constant and $B = B(x,y,z)$ is an arbitrary function. What I would really like to do is write a function which evaluates the above at some ordinate given by a label $rm{f}$, and another, $rm{s}$, and compute the difference between them. In other words: I would like Mathematica to analytically compute

$$ left[dfrac{1}{mu} dfrac{partial B}{partial x} right]^{rm f}_{rm s}:=dfrac{1}{mu_{rm f}}dfrac{partial B_{rm f}}{partial x}-dfrac{1}{mu_{rm s}}dfrac{partial B_{rm s}}{partial x}$$
particularly so that the variables with subscripts can be used later on (e.g. replaced with specific values, or manipulated as a function so that it can be obtained as a solution to a differential equation). In particular, I would like (a) the subscripting procedure to work inside of differential operators as seen above, and (b) the option to specify variables which are unchanged by the above procedure (e.g. due to being continuous across the interface).

I know you can use something like

AtFunc[func_[vars__], label_] := Subscript[func, label][vars]

to give a label to a function, but how can this be adapted to (a) apply to constants as well as functions, and (b) adjusted so that the labelling goes inside the derivatives?

EDIT: I am only looking for a solution which gives me a symbolic expression, not a numerical solution to any problem. The goal of this is to automate the "writing down" of the model I am to solve, before applying numerical (or asymptotic) methods of any sort.

diff[a_,b_]= B'[a] - B'[b]

diff[fun:B_,a_,b_]

 diff[Sin, 0, 1]
(*1 - Cos[1]*)