Thanks Ron and Sebastian for your help.  The code to calculate the RHS of the ODEs is really large so I can&#39;t really post it.  I am solving a reaction-diffusion-advection system in 2D.  The velocity (which contributes to the advection) is a complicated function of the concentrations of the species in the reaction-diffusion system.  At each time step, I have to solve a separate, time-independent ODE for the velocity, which is then incorporated into the RHS.  Because this ODE happens to be linear, I can solve it using a spectral method (for my boundary conditions of interest, regardless of my differencing scheme in the reaction-diffusion-advection system).  The presence of this velocity, manifest in the advection term, makes the Jacobian dense, even if I use low-order finite difference for my differencing.  Since I necessarily have a dense Jacobian, I typically use pseudospectral methods for computing the derivatives to get high accuracy with fewer grid points, though I have also implemented finite difference methods.<br>
<br>I have pasted below the first few lines of output of a profile I did using the magic function <span style="font-family: courier new,monospace;">%prun</span> with <span style="font-family: courier new,monospace;">ipython -pylab</span>.<br>
<br><span style="font-family: courier new,monospace;">      36447356 function calls in 1212.986 CPU seconds</span><br style="font-family: courier new,monospace;"><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">   Ordered by: internal time</span><br style="font-family: courier new,monospace;">
<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">   ncalls  tottime  percall  cumtime  percall filename:lineno(function)</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">        1  903.674  903.674 1212.976 1212.976 {scipy.integrate._odepack.odeint}</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">   678300  157.372    0.000  239.023    0.000 basic.py:358(_raw_fftnd)</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;">  1356600   24.225    0.000   35.504    0.000 basic.py:109(_fix_shape)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"></span><br>
<font face="arial,helvetica,sans-serif">So, aside from <span style="font-family: courier new,monospace;">odeint</span> itself, most of the computational effort is spent computing derivatives and solving for the velocity, all of which is just using 2D FFTs.  Curiously, <span style="font-family: courier new,monospace;">odeint</span> spends a lot of time doing something that does not require evaluation of the RHS of the ODEs.</font>  As I watch the integration proceed, it takes several time steps very rapidly, and then pauses at a particular time point, presumably to adjust the step size.  My guess is that Sebastian is right: this is slow because the Jacobian is large and dense.  Do you think just doing an explicit RK4 time stepping or something might help?  Or are there tricks we can play to get <span style="font-family: courier new,monospace;">odeint</span> to roll along more quickly?  I am passing an object with all sorts of attributes as an argument to the RHS function, though only some are used.  I am assuming the RHS function is all executed in Python and then the return values are all that undergo type conversion, right?  I&#39;m guessing it&#39;s not so much type conversion that&#39;s the big slowdown, or do you think I&#39;m mistaken?  <span style="font-family: courier new,monospace;">odeint<font face="arial,helvetica,sans-serif"> has been a nice little workhorse in my work so far!<br>
</font></span>