[R] solving ODE's in matrix form with lsoda()

[Jorge Ahumada]

Hello there,

Suppose you want to solve the following system of ODE's (a simple  
Lotka-Volterra predator prey model)

dP/dt = beta*P*V - mu*P
dV/dt = r*V - beta*P*V

where P and V are the numbers of predators and prey. Now, this is  
easy to do, but suppose you have a system of equations like this,

dP1/dt = beta1*P1*V1 - mu1*P1
dP2/dt = beta2*P2*V2 - mu2*P2
dV1/dt = r1*V1 - beta1*P1*V1
dV2/dt = r1*V2 - beta1*P2*V2

System 1 and system 2 are independent but that doesn't need to be the  
case. Now can you specify this system in lsoda() as:

dP/dt = beta*P*V - mu*P
dV/dt = r*V - beta*P*V

but now the initial state variables are 1x2 vectors and the  
parameters beta, mu and r are also vectors of size 1x2:

y=c(10,20,10,20)
parms = matrix(c(0.05,0.1,0.2,0.05,0.1,0.2),nc=3,byrow=T)

model = function(times,y,parms) {
y=matrix(y,nc=2,byrow=T)
beta=parms[,1]; mu = parms[,2], r = parms[,3]

dP/dt = beta*P*V - mu*P
dV/dt = r*V - beta*P*V

list(c(dP/dt,dV/dt)

}

I tried this but it does not give me the right answer. lsoda() gives  
warnings and produces weird results..

thanks, Jorge