We present BFComp, an automated framework based on Sum-Of-Squares (SOS) optimization and δ-decidability over the reals, to compute Bisimulation Functions (BFs) that characterize Input-to-Output Stability (IOS) of dynamical systems. BFs are Lyapunov-like functions that decay along the trajectories of a given pair of systems, and can be used to establish the stability of the outputs with respect to bounded input deviations.In addition to establishing IOS, BFComp is designed to provide tight bounds on the squared output errors between systems whenever possible. For this purpose, two SOS optimization formulations are employed: SOSP 1, which enforces the decay requirements on a discretized grid over the input space, and SOSP 2, which covers the input space exhaustively. SOSP 2 is attempted first, and if the resulting error bounds are not satisfactory, SOSP 1 is used to compute a Candidate BF (CBF). The decay requirement for the BFs is then encoded as a δ-decidable formula and validated over a level set of the CBF using the dReal tool. If dReal produces a counterexample containing the states and inputs where the decay requirement is violated, this pair of vectors is used to refine the input-space grid and SOSP 1 is iterated.By computing BFs that appeal to a small-gain theorem, the BFComp framework can be used to show that a subsystem of a feedback-composed system can be replaced-with bounded error-by an approximately equivalent abstraction, thereby enabling approximate model-order reduction of dynamical systems. We illustrate the utility of BFComp on a canonical cardiac-cell model, showing that the four-variable Markovian model for the slowly activating Potassium current IKs can be safely replaced by a one-variable HodgkinHuxley-type approximation.