chemo_1=function(t,u,parms){ # # Function chemo_1 computes the t derivative vectors of # u1(x,t), u2(x,t) # # One vector to two vectors u1=rep(0,nx);u2=rep(0,nx); for(i in 1:nx){ u1[i]=u[i]; u2[i]=u[i+nx]; } # # u1x, u2x u1x=rep(0,nx);u2x=rep(0,nx); u1x=dss004(xl,xu,nx,u1); u2x=dss004(xl,xu,nx,u2); # # BCs u1x[1]=0; u1x[nx]=0; u2x[1]=0; u2x[nx]=0; # # Nonlinear term u1u2x=rep(0,nx); for(i in 1:nx){ u1u2x[i]=2*u2[i]/u1[i]*u1x[i]; } # # u1u2xx, u2xx u2xx=rep(0,nx);u1u2xx=rep(0,nx); u2xx =dss004(xl,xu,nx, u2x); u1u2xx=dss004(xl,xu,nx,u1u2x); # # PDEs u1t=rep(0,nx);u2t=rep(0,nx); for(i in 1:nx){ u1t[i]=-k*u2[i]; u2t[i]=D*(u2xx[i]-u1u2xx[i]); } # # Two vectors to one vector ut=rep(0,2*nx); for(i in 1:nx){ ut[i] =u1t[i]; ut[i+nx]=u2t[i]; } # # Increment calls to chemo_1 ncall <<- ncall+1; # # Return derivative vector return(list(c(ut))); }