🚀 go-pugleaf

  As a appetizer to a more complex problem, I've been given the
task to efficiently solve the following infinite set of linear
ODEs (which arise when using an stochastic model to describe the
time behaviour of disease in a population of susceptibles and 
infecteds):

               _                                        _
dP(i,j,t)     | (i+1)(j-1)                  ij           |
--------- = A |----------- P(i+1,j-1,t) - ------ P(i,j,t)| 
dt            |_  i+j-1                    i+j-1        _|
               _                    _
              |                      |
          + B | P(i-1,j,t) - P(i,j,t)|
              |_                    _|
               _                            _
              |                              |
          + C | (j+1) P(i,j+1,t) - j P(i,j,t)|
              |_                            _|

               _                            _
              |                              |
          + D | (i+1) P(i+1,j,t) - i P(i,j,t)|
              |_                            _|


where  P(i,j,t) is the probability that the population has
i susceptibles and j infecteds at time t. A, B, C and D are
constants. Also P(i,j,t)=0 if i < 0 or j < 0.
i and j are integers, varying from 0 to infinity. To solve the
problem numerically we actually set the maximum values of 
i and j to imax=350 and jmax=150, respectively.
The initial conditions are:

P(N,1,0) = 1 , and P(i,j,0) = 0 for other values of i and j.

N is the initial size of the population (N=100).

So far we have solved this problem using the routine ODEINT
from "Numerical Recipes in C". It takes about 15 hours on
an Sparcstation (time range from 0 to 300 in steps of 5).


Currently I'm coding the problem above to use the fortran
routines available in the package VODE. 
Notice that we have (imax + 1)*(jmax + 1) linear ODEs.

I wondered if someone has any suggestions to tackle a problem
like this, for example, by rewritting the ODEs in a different
form, by using another package (fortran or C), by using Maple
or Mathematica, etc.

Thanks for any suggestions.

   Roque
   oliveria@engin.umich.edu