diff --git a/R/setup.R b/R/setup.R
index 0a9076ad..c3589f2f 100644
--- a/R/setup.R
+++ b/R/setup.R
@@ -78,6 +78,8 @@ xde_setup = function(modelName,
   pars$compute = list()
   class(pars$compute) = "xde"
 
+  pars$outputs = list()
+
   return(pars)
 }
 
@@ -116,7 +118,6 @@ xde_setup_mosy = function(modelName,
 
                      # forcing
                      kappa=NULL
-
 ){
 
   pars = make_parameters_xde()
@@ -142,6 +143,8 @@ xde_setup_mosy = function(modelName,
   pars$compute = list()
   class(pars$compute) = "na"
 
+  pars$outputs = list()
+
   return(pars)
 }
 
@@ -177,6 +180,8 @@ xde_setup_aquatic = function(modelName,
   pars$compute = list()
   class(pars$compute) = "na"
 
+  pars$outputs = list()
+
   return(pars)
 }
 
@@ -236,6 +241,8 @@ xde_setup_human = function(modelName,
   pars$compute = list()
   class(pars$compute) = "human"
 
+  pars$outputs = list()
+
   return(pars)
 }
 
@@ -282,5 +289,7 @@ xde_setup_cohort = function(modelName, F_eir,
   pars$compute = list()
   class(pars$compute) = "cohort"
 
+  pars$outputs = list()
+
   return(pars)
 }
diff --git a/R/solving.R b/R/solving.R
index f474e95d..9d77107f 100644
--- a/R/solving.R
+++ b/R/solving.R
@@ -23,7 +23,7 @@ xde_stable_orbit = function(pars, Ymax=10){
   steady$terms = pars$orbits
   for(i in 1:length(steady$terms))
     steady$terms[[i]] = tail(steady$terms[[i]],365)
-  pars$stable_orbits <- steady
+  pars$outputs$stable_orbits <- steady
   return(pars)
 }
 
@@ -36,7 +36,7 @@ xde_steady = function(pars){
   y0 = get_inits(pars)
   pars1 <- dde2ode_MYZ(pars)
   rootSolve::steady(y=y0, func = xDE_diffeqn, parms = pars1)$y -> y_eq
-  pars$steady = parse_deout_vec(y_eq, pars)
+  pars$outputs$steady = parse_deout_vec(y_eq, pars)
   return(pars)
 }
 
@@ -49,7 +49,7 @@ xde_solve.ode = function(pars, Tmax=365, dt=1){
   tt = seq(0, Tmax, by=dt)
   y0 = get_inits(pars)
   deSolve::ode(y = y0, times = tt, func = xDE_diffeqn, parms = pars, method = "lsoda") -> out
-  pars$orbits = parse_deout(out, pars)
+  pars$outputs$orbits = parse_deout(out, pars)
   return(pars)
 }
 
@@ -62,7 +62,7 @@ xde_solve.dde = function(pars, Tmax=365, dt=1){
   tt = seq(0, Tmax, by=dt)
   y0 = get_inits(pars)
   deSolve::dede(y = y0, times = tt, func = xDE_diffeqn, parms = pars, method = "lsoda") -> out
-  pars$orbits = parse_deout(out, pars)
+  pars$outputs$orbits = parse_deout(out, pars)
   return(pars)
 }
 
@@ -75,7 +75,7 @@ xde_solve.aqua = function(pars, Tmax=365, dt=1){
   tt = seq(0, Tmax, by=dt)
   y0 = get_inits_L(pars)
   deSolve::ode(y = y0, times = tt, func = xDE_diffeqn_aquatic, parms = pars, method = "lsoda") -> out
-  pars$orbits = parse_deout(out, pars)
+  pars$outputs$orbits = parse_deout(out, pars)
   return(pars)
 }
 
@@ -88,7 +88,7 @@ xde_solve.aqua_dde = function(pars, Tmax=365, dt=1){
   tt = seq(0, Tmax, by=dt)
   y0 = get_inits_L(pars)
   deSolve::dede(y = y0, times = tt, func = xDE_diffeqn_aquatic, parms = pars, method = "lsoda") -> out
-  pars$orbits = parse_deout(out, pars)
+  pars$outputs$orbits = parse_deout(out, pars)
   return(pars)
 }
 
@@ -101,7 +101,7 @@ xde_solve.mosy = function(pars, Tmax=365, dt=1){
   tt = seq(0, Tmax, by=dt)
   y0 = get_inits(pars)
   deSolve::ode(y = y0, times = tt, func = xDE_diffeqn_mosy, parms = pars, method = "lsoda") -> out
-  pars$orbits = parse_deout(out, pars)
+  pars$outputs$orbits = parse_deout(out, pars)
   return(pars)
 }
 
@@ -114,7 +114,7 @@ xde_solve.mosy_dde = function(pars, Tmax=365, dt=1){
   tt = seq(0, Tmax, by=dt)
   y0 = get_inits(pars)
   deSolve::dede(y = y0, times = tt, func = xDE_diffeqn_mosy, parms = pars, method = "lsoda") -> out
-  pars$orbits = parse_deout(out, pars)
+  pars$outputs$orbits = parse_deout(out, pars)
   return(pars)
 }
 
@@ -127,7 +127,7 @@ xde_solve.human = function(pars, Tmax=365, dt=1){
   tt = seq(0, Tmax, by=dt)
   y0 = get_inits(pars)
   deSolve::ode(y = y0, times = tt, func = xDE_diffeqn_human, parms = pars, method = "lsoda") -> out
-  pars$orbits = parse_deout(out, pars)
+  pars$outputs$orbits = parse_deout(out, pars)
   return(pars)
 }
 
@@ -141,6 +141,6 @@ xde_solve.cohort = function(pars, Tmax=365, dt=1){
   tt = seq(0, Tmax, by=dt)
   y0 = get_inits(pars)
   deSolve::ode(y = y0, times = tt, func = xDE_diffeqn_cohort, parms=pars, F_eir = pars$F_eir, method = "lsoda") -> out
-  pars$orbits = parse_deout(out, pars)
+  pars$outputs$orbits = parse_deout(out, pars)
   return(pars)
 }
diff --git a/vignettes/adult_RM.Rmd b/vignettes/adult_RM.Rmd
index 610850d6..b8e712be 100644
--- a/vignettes/adult_RM.Rmd
+++ b/vignettes/adult_RM.Rmd
@@ -277,7 +277,7 @@ Now, we can solve the equations using `xde_solve` and compare the output to what
 
 ```{r xde_solve}
 xde_solve(MYZeg, Tmax=50, dt=1) -> MYZeg 
-out2 <- MYZeg$orbits$deout
+out2 <- MYZeg$outputs$orbits$deout
 sum(abs(out1-out2))==0 
 ```
 
diff --git a/vignettes/aqua_basic.Rmd b/vignettes/aqua_basic.Rmd
index 339d2273..0ad07c28 100644
--- a/vignettes/aqua_basic.Rmd
+++ b/vignettes/aqua_basic.Rmd
@@ -139,7 +139,7 @@ Mo = list(Gm = c(250, 500, 170))
 
 ```{r}
 xde_setup_aquatic("aqbasic", nHabitats = 3, Lname = "basic", Lopts = Lo, MYZopts = Mo) -> aqbasic
-xde_solve(aqbasic, Tmax=50, dt=1)$deout -> out2
+xde_solve(aqbasic, Tmax=50, dt=1)$output$orbits$deout -> out2
 sum(abs(out1-out2)) == 0
 ``` 
 
diff --git a/vignettes/aqua_trace.Rmd b/vignettes/aqua_trace.Rmd
index abcf529e..fa427b77 100644
--- a/vignettes/aqua_trace.Rmd
+++ b/vignettes/aqua_trace.Rmd
@@ -168,7 +168,7 @@ xde_setup_mosy("mosy1",
 ```
 
 ```{r}
-xde_solve(mosy1,Tmax=50,dt=1)$orbits$deout -> out2
+xde_solve(mosy1,Tmax=50,dt=1)$outputs$orbits$deout -> out2
 sum(out2-out1)
 ```
 
diff --git a/vignettes/ex_534.Rmd b/vignettes/ex_534.Rmd
index 36e281ea..e2eb61e5 100644
--- a/vignettes/ex_534.Rmd
+++ b/vignettes/ex_534.Rmd
@@ -337,7 +337,7 @@ We can test the idea:
 ```{r}
 mod534$calU = params$calU
 mod534 <- xde_solve(mod534,Tmax=365,dt=1)
-mod534$orbits$deout -> out2
+mod534$outputs$orbits$deout -> out2
 sum(abs(out2-out1))
 ```
 
diff --git a/vignettes/human_sip.Rmd b/vignettes/human_sip.Rmd
index 96955608..11ef0d75 100644
--- a/vignettes/human_sip.Rmd
+++ b/vignettes/human_sip.Rmd
@@ -139,7 +139,7 @@ xde_setup_cohort("test_SIP", F_eir1, "SIP", HPop=Hpop, Xopts = Xo) -> test_SIP
 ```
 
 ```{r}
-xde_solve(test_SIP, 365, 365)$deout -> out2
+xde_solve(test_SIP, 365, 365)$outputs$orbits$deout -> out2
 approx_equal(out2,out1) 
 ```
 
diff --git a/vignettes/human_sis.Rmd b/vignettes/human_sis.Rmd
index a424d1e5..964d84d5 100644
--- a/vignettes/human_sis.Rmd
+++ b/vignettes/human_sis.Rmd
@@ -138,6 +138,6 @@ xde_setup_cohort("test_SIS", F_eir1, "SIS", HPop=Hpop, Xopts = Xo) -> test_SIS
 ```
 
 ```{r}
-xde_solve(test_SIS, 365, 365)$deout -> out2
+xde_solve(test_SIS, 365, 365)$outputs$orbits$deout -> out2
 approx_equal(out2,out1)
 ```
diff --git a/vignettes/vc_lemenach.Rmd b/vignettes/vc_lemenach.Rmd
index 7743519e..aeb76830 100644
--- a/vignettes/vc_lemenach.Rmd
+++ b/vignettes/vc_lemenach.Rmd
@@ -170,7 +170,7 @@ itn_mod = setup_itn_lemenach(itn_mod, phi = ITN_cov)
 
 ```{r}
 xde_solve(itn_mod, 5*365) -> itn_mod
-itn_mod$orbits$deout -> out
+itn_mod$outputs$orbits$deout -> out
 ```
 
 ## Plot