Dear David, It probably would be easier to pass a function that, given the current states and the parameters, returned a Float64 vector of rates. The main issue that I have is that eval can't see all the variables defined in the function scope, so my code runs fine if run globally, but not as a function. If I go along that route, it would probably be best to define an appropriate model type, rather like Tim Vaughan's MASTER program:
https://github.com/CompEvol/MASTER/wiki/Trajectory-specification I'll tweak my code to work with functions, and push to github when I get a chance. Best Simon