In financial and egineering problems, we are often faced with solving Partial-Integro Differential Equations (PIDEs). Rarely we can find an analytic solution in a closed form expression for these PIDEs, hence we turn to numerical schemes to accurately approximate the solution instead. Classically these methods are based on finite difference methods, however, we can turn certain kinds of PIDEs into a probabilistic representation, called Forward Backward Stochastic Differential Equations with Jumps (FBSDEJs). Solving the PIDE can now be done alternatively by solving a FBSDEJ. In this thesis we will first investigate the stochastic framework behind FBSDEJs and we will look into the uniqueness and existence of their solutions. Furthermore we propose a new numerical method which can efficiently solve FBSDEJs. The semi-discretisation is based on the classical Backward Differentiation Formula (BDF) methods, for the computation of the conditional expectations we use the COS method which makes use of Fourier cosine expansions, exploiting the knowledge we have about characteristic functions. Finally we implement the new method and we investigate it extensively both numerically and theoretically. We show that the BDFn schemes are highly stable and efficient for computing FBSDEJs, the initial steps still have to be investigated in greater detail so that we can make use of the high-order BDFn schemes.