times <- seq(-500,500,5)
Ber78Obl <- sapply(times, function(t) obliquity(t, ber78))
La04Obl <- sapply(times, function(t) obliquity(t, la04))
plot(times, Ber78Obl)
lines(times, La04Obl)
plot(times, (Ber78Obl - La04Obl)*180/pi)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.