rkf45=function(nt,h,t,y) { # # Function rkf45 implements a fourth order Runge Kutta method # embedded in a fifth order Runge Kutta method. The ODE routine # has to be renamed for a new application. # # The arguments are # # Input # # nt - number of rkf45 steps between the starting and final # points along the solution # # h - integration step (constant) # # t - initial value of the independent variable # # y - initial dependent variable vector # # Output # # y - solution vector after nt integration steps # # nt rkf45 steps for(i2 in 1:nt){ if(nint==1){ yb=y; tb=t; rk1=apoptosis_2(tb,yb)*h y=yb+0.25*rk1; t=tb+0.25*h; rk2=apoptosis_2(t,y)*h y=yb+(3/32)*rk1+(9/32)*rk2; t=tb+(3/8)*h; rk3=apoptosis_2(t,y)*h y=yb+(1932/2197)*rk1-(7200/2197)*rk2+(7296/2197)*rk3; t=tb+(12/13)*h; rk4=apoptosis_2(t,y)*h y=yb+(439/216)*rk1-8*rk2+(3680/ 513)*rk3-( 845/4104)*rk4; t=tb+h; rk5=apoptosis_2(t,y)*h y=yb-(8/27)*rk1+2*rk2-(3544/2565)*rk3+(1859/4104)*rk4-(11/40)*rk5; t=tb+0.5*h; rk6=apoptosis_2(t,y)*h y=yb+(16/135)*rk1+(6656/12825)*rk3+(28561/56430)*rk4-(9/50)*rk5+ (2/55)*rk6; t=tb+h; } if(nint==2){ yb=y; tb=t; rk1=apoptosis_2(tb,yb)*h y=yb+0.25*rk1; t=tb+0.25*h; rk2=apoptosis_2(t,y)*h y=yb+(3/32)*rk1+(9/32)*rk2; t=tb+(3/8)*h; rk3=apoptosis_2(t,y)*h y=yb+(1932/2197)*rk1-(7200/2197)*rk2+(7296/2197)*rk3; t=tb+(12/13)*h; rk4=apoptosis_2(t,y)*h y=yb+(439/216)*rk1-8*rk2+(3680/ 513)*rk3-( 845/4104)*rk4; t=tb+h; rk5=apoptosis_2(t,y)*h y=yb-(8/27)*rk1+2*rk2-(3544/2565)*rk3+(1859/4104)*rk4-(11/40)*rk5; t=tb+0.5*h; rk6=apoptosis_2(t,y)*h # # Fourth order step y4=yb+(25/216)*rk1+( 1408/2565)*rk3+( 2197/4104)*rk4-( 1/5)*rk5; # # Fifth order step y=yb+(16/135)*rk1+(6656/12825)*rk3+(28561/56430)*rk4-(9/50)*rk5+ (2/55)*rk6; t=tb+h; # # Truncation error estimate ee=y-y4 ee <<- ee } } return(c(y)) }