🚀 go-pugleaf

I am trying to solve a system of simultaneous PDE's and ODE's.  

There are 3 PDE's: standard one dimensional advection-diffusion equations with a nonlinear reaction (uptake) term.  The reaction term couples the PDE's to a system (at each spatial node) of 28 ODE's.  The ODE's come from a biological growth model and are
thus nonlinear and very complicated.

I would welcome any and all suggestions (with references to the literature, if 
available) as to numerical algorithms to solve these equations.

One approach I am considering is to apply a finite difference approximation to the right hand sides of the PDE's, thus converting them to ODE's and solving the
resultant system of 31n (n = # of nodes) ODE's with an adaptive stepsize method
of some sort. Comments?

I have some familiarity with basic numerical methods.  My preference would be to 
use 'canned' FORTRAN numerical routines (ala netlib) to the maximum extent possible.
My hardware choices include PC (486-33), CRAY, SUN (SPARC10), in order of familiarity and ease of use.

Thanks in advance for all the useful help I am sure I will get. Responses can be
posted here or E-mailed directly to me, whichever seems more appropriate.