drug_2=function(t,U,parms){ # # Function drug_2 computes the t derivative vector # of the u,v,s vectors # # One vector to three vectors u=rep(0,nx);v=rep(0,nx); s=rep(0,nx); for(i in 1:nx){ u[i]=U[i]; v[i]=U[i+nx]; s[i]=U[i+2*nx]; } # # ux, sx # ux=dss004(xl,xu,nx,u); # sx=dss004(xl,xu,nx,s); # # Boundary conditions ux=rep(0,nx);sx=rep(0,nx); ux[1]=-(kr/D)*(ua-u[1]); ux[nx]=(kr/D)*(ua-u[nx]); sx[1]=0;sx[nx]=0; # # uxx, sxx nl=2;nu=2; uxx=dss044(xl,xu,nx,u,ux,nl,nu); sxx=dss044(xl,xu,nx,s,sx,nl,nu); # # PDEs ut=rep(0,nx);vt=rep(0,nx) st=rep(0,nx); for(i in 1:nx){ ut[i]=D*uxx[i]+E*sxx[i]+f_u(ub,vb,u[i],v[i]); vt[i]=g_v(ub,vb,u[i],v[i]); st[i]=alpha*u[i]-beta*s[i]+gamma*ut[i]; } # # Three vectors to one vector Ut=rep(0,3*nx); for(i in 1:nx){ Ut[i] =ut[i]; Ut[i+nx] =vt[i]; Ut[i+2*nx]=st[i]; } # # Increment calls to drug_2 ncall <<- ncall+1; # # Return derivative vector return(list(c(Ut))); }