Halflife for decay function

Schrick-Noah_Homework-3.R

@@ -140,10 +140,35 @@ y.init <- 10000
 
 decay.ode45.sol <- ode45(f=function(t,y){decay.f(t,y,k=k)},
                          y=y.init, t0=tmin, tfinal=tmax)
-plot(decay.ode45.sol$t,decay.ode45.sol$y)
+plot(decay.ode45.sol$t,decay.ode45.sol$y, xlab="t", ylab="y",
+     main="Numerical Decay Model Solution Using Pracma ODE45")
 par(new=T)
 lines(seq(tmin,tmax,len=50),
-      y.init*exp(-k*seq(tmin,tmax,len=50)))
+      y.init*exp(-k*seq(tmin,tmax,len=50)), col="red")
+abline(h=y.init/2, lty="dashed", lwd=3)
+
+#Solve for half-life with bisection
+findZeroBisect <- function(func, xl, xr, tol, maxsteps=1e6){
+  steps <- 0
+  if (func(xl) * func(xr) >= 0){
+    stop("Bad interval: try again with different endpoints")
+  }
+  
+  repeat{
+    xm <- (xl + xr)/2
+    if(abs(func(xm)) < tol || steps >= maxsteps){
+      break
+    }
+    
+    ifelse(func(xl) * func(xm) > 0, xl<-xm, xr<-xm)
+    steps <- steps + 1
+  }
+  return(c(xm, func(xm), steps))
+}
+fun.hl <- function(x) {-5000+(y.init)*exp(-k*x)}
+
+halflife <- findZeroBisect(fun.hl, 5000, 6500, 1e-8)[1]
+abline(v=halflife, lty="dotted", lwd=2, col="blue")
 
 ## 4. Solve the predator-prey model numerically