Interesting. Can you get a simplified example which recreates this problem and post it in an issue to DiffEqJump.jl? https://github.com/JuliaDiffEq/DiffEqJump.jl/issues . It looks like the first jump is re-using the same random number each time. We have tests for that when doing looping but might have omitted the test when using a MonteCarloProblem, but that is obviously a good think to fix here.

If you don’t mind, at first I post the code here. Maybe I do something wrong when I convert the solution into array.

using DifferentialEquations, Plots
rate(u,p,t) = 2
affect!(integrator) = (integrator.u[1] = integrator.u[1]/2)
jump = VariableRateJump(rate,affect!)
jump_prob = JumpProblem(prob,Direct(),jump)
monte_prob = MonteCarloProblem(jump_prob)
sol = solve(monte_prob,SRIW1(),num_monte=3,parallel_type=:none,dt=0.001,adaptive=false)
#I need the trajectories of equal lengths, so I do that
A = Array{Float64}(undef, 1000, 3)
for i in 1:3
A[:,i] = Array(sol[i])[1,1:1000]
end
pyplot()
p = plot(A, legend = false, title="$proc_name $n")
display(p)

Oh, the jumps add a time point in there themselves since they introduce a discontinuity. If you don’t want that save, you do save_positions=(false,false) in VariableRateJump. Then Array(sol) will work and line up with the time points you’re expecting. I am not sure if the first jump is caused actually the same or if this rescaling you did made them look the same, so give this a shot first.

Yup okay, this requires an issue. Somehow the callback init which is supposed to re-sample that random number isn’t being called. I think it’s confined to only variable rate jumps though.

What the sense in this line for the jumps? affect!(integrator) = (integrator.u[1] = 1) or integrator.u[1] = integrator.u[1]
As I understand, it affects on the size of jumps at time points, doesn’t it? But in what time points?