diff --git a/Schrick-Noah_Homework-4.R b/Schrick-Noah_Homework-4.R
index 7085c18..d41e96b 100644
--- a/Schrick-Noah_Homework-4.R
+++ b/Schrick-Noah_Homework-4.R
@@ -259,11 +259,6 @@ plot(best.sol$x, best.sol$y[,1], type="l")
 
 ## 4. Position of the earth and moon
 # a. Plotly
-G <- 6.673e-11   # m^3 kg-1 s^-2
-M_S <- 1.9891e30 # sun kg
-M_E <- 5.98e24   # earth kg
-M_m <- 7.32e22   # moon kg
-mu_sun <- G*M_S*(86400^2)/1e9 #  km^3/days^2
 
 # t0: jan 1, 1999 00:00:00am
 x0_earth <- c(-27115219762.4, 132888652547.0, 57651255508.0)/1e3 # km
@@ -317,15 +312,39 @@ mean(sqrt(rowSums(cbind(sunearth.sol$y[,1]^2,
                         sunearth.sol$y[,3]^2)))) 
 # 150 million km
 
-mean(sqrt(rowSums(cbind((sunearth.sol$y[,1] - sunmoon.sol$y[,1])^2, 
-                        (sunearth.sol$y[,2] - sunmoon.sol$y[,2])^2, 
-                        (sunearth.sol$y[,3] - sunmoon.sol$y[,3])^2))))
+sunmoon.f <- function(t,y){
+  G <- 6.673e-11 # m^3 kg^-1 s^-2
+  M_E <- 5.98e24
+  M_m <- 7.32e22  
+  mu_sun <- G*M_E*(86400^2)/1e9  # km^3/days^2
+  
+  x_moon <- y[1:3]  # xyz position of moon wrt earth
+  v_moon <- y[4:6]  # xyz velocity of moon, also part of ydot
+  r_moon <- sqrt(sum(x_moon^2))  # distance moon to earth
+  acc_moon <- -mu_sun*x_moon/r_moon^3  # part of ydot
+  
+  ydot <- c(v_moon, acc_moon)
+  return(matrix(ydot))
+}
 
-mean(sqrt(rowSums(cbind(sunmoon.sol$y[,1]^2, 
-                        sunmoon.sol$y[,2]^2, 
-                        sunmoon.sol$y[,3]^2))))
-# Not correct: giving 60m km - should be ~380k
-# Why: Don't think M_m is used in function
+x0_moonearth <- x0_moon - x0_earth
+v0_moonearth <- v0_moon - v0_earth
+y.init.moonearth <- c(x0_moonearth, v0_moonearth)
+earthmoon.sol <- rk4sys(f=sunmoon.f, a=0,b=365.25, y0=y.init.moonearth, n=365)
+earthmoon.df <- data.frame(x=earthmoon.sol$y[,1],
+                         y=earthmoon.sol$y[,2],
+                         z=earthmoon.sol$y[,3])
+
+mean(sqrt(rowSums(cbind(earthmoon.sol$y[,1]^2, 
+                        earthmoon.sol$y[,2]^2, 
+                        earthmoon.sol$y[,3]^2))))
+
+fig <- plot_ly(earthmoon.df, x = ~x, y = ~y, z = ~z, name="moon", type = "scatter3d", mode="markers") 
+fig <- fig %>% layout(scene = list(xaxis = list(title = 'x'),
+                                   yaxis = list(title = 'y'),
+                                   zaxis = list(title = 'z')))
+fig <- add_trace(fig, x = 0, y = 0, z = 0, mode="markers", color = I("green"), name="earth")
+fig
 
 # b. Find eclipses
 norm_vec <- function(x) sqrt(rowSums(cbind((x^2))))