# Make sure that your first graph is still open, # and that the original graph was drawn so that # the y-axis goes from 0 to 3.5 y_vals = c(0, 3.5) # Vector of initial conditions y_0 = seq(0.2, 2, by = 0.2) # For-loop for (i in 1:length(y0)) { # Re-set start each time start = c(x_0, y_0[i]) # Solve the ODE output = as.data.frame(rk4(start, times, odef, parms)) # Plot our values on the same graph par(new = TRUE) plot(output$"2", type = "l", col = c(i), ylab = "", xlab = "", ylim = y_vals, xaxt = "n", yaxt = "n") }