bz_1=function(t,u,parms){ # # Function bz_1 computes the t derivative vector # of the u1,u2 vectors # # 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=dss004(xl,xu,nx,u1); u2x=dss004(xl,xu,nx,u2); # # Boundary conditions u1x[1]=0;u1x[nx]=0; u2x[1]=0;u2x[nx]=0; nl=2;nu=2; # # u1xx, u2xx u1xx=dss044(xl,xu,nx,u1,u1x,nl,nu); u2xx=dss044(xl,xu,nx,u2,u2x,nl,nu); # # PDEs u1t=rep(0,nx);u2t=rep(0,nx); for(i in 1:nx){ u1t[i]=u1xx[i]+u1[i]*(1-u1[i]-r*u2[i]); u2t[i]=u2xx[i]-b*u1[i]*u2[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 bz_1 ncall <<- ncall+1; # # Return derivative vector return(list(c(ut))); }