I think an additional problem is that the derivatives function must be 'told' where to 'look' in y for the values of R, C, and P that are used to compute Rprime,Cprime, Pprime - it has no way of 'knowing' that the model's state vector y consists of the variables (R,C,P). Appended below is an example of code that works to solve a similar system of equations. Perhaps it might be useful to add something like that to the ?lsoda page?
> function(times,y,p) > { > Rprime <- > (R*(1-R))-((xc*yc*C*R)/(R+R0))-((w*xp*ypr*P*R)/(R02+((1-w)*C)+(w*R))) > Cprime <- > (-1*(xc*C)*(1-(yc*R)/(R+R0)))-(((1-w)*xp*ypc*P*C)/((w*R)+((1-w)*C)+C0)) > Pprime <- > >(-1*P)-(((1-w)*xp*ypc*C*P)/((w*R)+((1-w)*C)+C0))+((w*xp*ypr*P*R)/((w*R)+((1-w)*C)+R02)) >list(c(Rprime, Cprime, Pprime)) >} require(odesolve); # Set parameter values bmaxc=3.3; bmaxb=2.25; Kc=4.3; Kb=15; ni=80; m=0.05; lambda=0.4; epsilon=0.25; delta=0.95; f2k=function(t,y,parms){ dy=rep(0,4); n=y[1]; c=y[2]; r=y[3]; b=y[4]; fcn=bmaxc*n/(Kc+n); fbc=bmaxb*c/(Kb+c); dy[1]=delta*(ni-n)-fcn*c; dy[2]=fcn*c-fbc*b/epsilon-delta*c; dy[3]=fbc*r-(delta+m+lambda)*r; dy[4]=fbc*r-(delta+m)*b; list(dy); } y0=c(1,70,.01,.01); times=(0:240)/4; parms=0; out=lsoda(y0, times, f2k, parms, rtol=1e-4, atol=1e-4); matplot(times,cbind(out[,3]/5,out[,5]),type="l",col=c("green","red"), lty=1, xlab="Time (d)", ylab="Prey (green) and Predators (red)"); Stephen P. Ellner ([EMAIL PROTECTED]) Department of Ecology and Evolutionary Biology Corson Hall, Cornell University, Ithaca NY 14853-2701 Phone (607) 254-4221 FAX (607) 255-8088 ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help