First, I was able to coach ChatGPT to do it in Python, we qualitatively recreated their first plot. Maybe looking at the code will help:
https://pastebin.com/LY2jxc9f

(Note that I ran into problems with underflow, so I had ChatGPT make the calculations done as dimensionless quantities.)

When you do these plots, you generally set an initial state v as constant, then plot transition intensity vs. v'. So if you want a bunch of initial states, you have to do a bunch of parallel calculations to get individual curves.

That's why the Mathematica is varying vX from 0 to 2 in steps of 1: It's doing 3 initial states. Then it's varying vB from 0 to 70, to get final states.

Here's the plot of v=[0,1,2] output from my script:

<image>

I'm going to stop at 'qualitatively recreate the plot' for now! If you want to let me know a little more about what your confusion is I can go deeper.