# Tutorial on Coupling Jumps

**URL:** https://discourse.julialang.org/t/tutorial-on-coupling-jumps/70957
**Category:** Modelling & Simulations
**Tags:** diffeq, sciml
**Created:** [November 4, 2021, 2:42pm UTC](https://discourse.julialang.org/t/tutorial-on-coupling-jumps/70957 "2021-11-04T14:42:36Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![tom-plaa](https://avatars.discourse-cdn.com/v4/letter/t/d26b3c/32.png) [@tom-plaa](https://discourse.julialang.org/u/tom-plaa)
#### Post date: [November 4, 2021, 2:42pm UTC](https://discourse.julialang.org/t/tutorial-on-coupling-jumps/70957/1 "2021-11-04T14:42:36Z")

</div>

Didn’t create an issue on the repo because this might be just my interpretation being wrong.  
I am trying to run code based on the tutorial at [Jump Diffusion page](https://diffeq.sciml.ai/stable/tutorials/jump_diffusion/#Coupling-Jump-Problems).

It mentions that this is supposed to be a Cox process (doubly stochastic) whose rate switches between 10 and 0, but I’m having trouble understanding how this code is supposed to achieve this.

The system that seems to be implied is something like  
x(t) = N\_a(y, t)  
y(t) = N\_c(t) - N\_b(t)  
z(t) = N\_b(t) - N\_c(t)

Where  
N\_a (y, t) ~ Poisson(10yt)  
N\_b(y, t) ~ Poisson(yt)  
N\_c(z, t) ~ Poisson(zt)

First of all, the code fails if I try to solve this system with almost any u0 except (1,1,1) or (1,0,0) with the error message

> “Tried to add a tstop that is behind the current time. This is strictly forbidden”

And here is what I get if I initialize to (1,1,1):

 ![Screenshot from 2021-11-04 14-41-53](https://global.discourse-cdn.com/julialang/original/3X/7/4/744d02bf2b3464ce1fd388b3d6d31df97fe3c3d8.png)

The code used, below:

```julia
rate(u,p,t) = u[2]*10
affect!(integrator) = integrator.u[1] += 1.
jump1 = ConstantRateJump(rate,affect!)
rate(u,p,t) = u[2]
affect!(integrator) = (integrator.u[2] -= 1.;integrator.u[3] += 1.)
jump2 = ConstantRateJump(rate,affect!)

rate(u,p,t) = u[3]
affect!(integrator) = (integrator.u[2] += 1.;integrator.u[3] -= 1.)
jump3 = ConstantRateJump(rate,affect!)
prob = DiscreteProblem(ones(1),(0.0,10.0))
jump_prob = JumpProblem(prob,Direct(),jump1,jump2,jump3)
plot(solve(jump_prob))

```