# Project 3 for the University of Tulsa's CS-7863 Sci-Stat Course
# Numerical Ordinary Differential Equations
# Professor: Dr. McKinney, Spring 2023
# Noah L. Schrick - 1492657

## 1. Approximate deriv. of sin(x) at x=pi from h=10^-1 -> 10^-20
forward.approx <- function(func, x, h){
  approx <- (func(x+h)-func(x))/h
}

forward.approx.table <- matrix(nrow = 0, ncol = 3)  
colnames(forward.approx.table) <- c("h", "approximation", "error")
for(h in 10^(seq(-1,-20,-1))){
  approx <- forward.approx(sin, pi, h)
  error <- abs(cos(pi)-approx)
  forward.approx.table <- rbind(forward.approx.table, c(h, approx, error)) 
}

plot(abs(log(forward.approx.table[,"h"],10)), forward.approx.table[,"error"],
     xlab="h [1e-x]", ylab="error (logscale)", type="o", log="y",
     main = "Forward Approximation of sin(x) at x=pi")


# Repeat with central difference approx
central.diff.approx <- function(func, x, h){
  approx <- (func(x+h)-func(x-h))/(2*h)
}
central.diff.table <- matrix(nrow = 0, ncol = 3)  
colnames(central.diff.table) <- c("h", "approximation", "error")
for(h in 10^(seq(-1,-20,-1))){
  approx <- central.diff.approx(sin, pi, h)
  error <- abs(cos(pi)-approx)
  central.diff.table <- rbind(central.diff.table, c(h, approx, error)) 
}

plot(abs(log(central.diff.table[,"h"],10)), central.diff.table[,"error"],
     xlab="h [1e-x]", ylab="error (logscale)", type="o", log="y",
     main = "Central Difference Approximation of sin(x) at x=pi")

## 2.
# a) Runge-Kutta

# b) Numerically solve the decay ode and compare to Euler and RK error

## 3. Use library function ode45 to solve the decay numerically

## 4. Solve the predator-prey model numerically

# a) k1=0.01, k2=0.1, k3=0.001, k4=0.05, prey(0)=50, pred(0)=15. t=0 -> 200

# plot

# b) Use Euler and RK to solve with h=10

# plot comparing Prey solutions to ode45

# c) Use k3=0.02

## 5. Solve SIR model numerically from t=0 -> 20

# a) a=0.5, b=1, S(0)=0.9, I(0)=0.1, R(0)=0

# plot

# b) a=3

## 6. Decomp of dinitrogen pentoxygen into nitrogen dioxide and molecular oxygen

# a) k1=1.0, k2=0.5, k3=0.2,k4=1.5, [N2O5]o=1, all other IC's=0, t=0 -> 10

# b) Increase k4 to make the intermediate species -> 0