wxMaxima(obj, out, tex = FALSE)

#components(determinant, Sigma inverse, argument for exponent) for joint normal distribution
eq_c <- c("Tw ~ ((((PH) + (tw)) * (ta - Tc + 2) + (1 + (tw)) * (Ec/w - 2/w) -(1 + (PH))) +
sqrt((((PH) + (tw)) * (ta - Tc + 2) + (1 +(tw)) * (Ec/w - 2/w) - (1 + (PH)))^2 - 4 * (1 + (PH) +
(tw)) *(-(PH) * (ta - Tc + 2) + (1 - (tw) * (ta - Tc + 2)) * (2/w -Ec/w))))/(2 * (1 + (PH) +
(tw)))",
"Tf1 ~ (th1) * (ta - (((((PH) + (tw)) * (ta - Tc + 2) + (1 + (tw)) *(Ec/w - 2/w) - (1 + (PH))) +
 sqrt((((PH) + (tw)) * (ta -Tc + 2) + (1 + (tw)) * (Ec/w - 2/w) - (1 + (PH)))^2 - 4 *(1 + (PH) +
  (tw)) * (-(PH) * (ta - Tc + 2) + (1 - (tw) *(ta - Tc + 2)) * (2/w - Ec/w))))/(2 * (1 + (PH) +
  (tw)))) -Tc + 2) - 1",
"Ef1 ~ (ph1)/(PH) * (w * (((((PH) + (tw)) * (ta - Tc + 2) + (1 +(tw)) * (Ec/w - 2/w) -
(1 + (PH))) + sqrt((((PH) + (tw)) *(ta - Tc + 2) + (1 + (tw)) * (Ec/w - 2/w) - (1 + (PH)))^2 -4 *
 (1 + (PH) + (tw)) * (-(PH) * (ta - Tc + 2) + (1 - (tw) *(ta - Tc + 2)) * (2/w - Ec/w))))/(2 *
 (1 + (PH) + (tw)))) -Ec + 2) - 1")
parl <- c("tw","PH","th1","ph1")
para_cont <- get_par(parl, eq_c)
cheqs0 <- para_cont$cheqs0
npar <- get_npar(cheqs0)
neq <- length(cheqs0)
sdv <- rep(NA, neq)
mv <- rep(NA, neq)
sigma <- expand.grid(1:neq, 1:neq)%>%apply(., 1, function(x)paste0(x, collapse=','))%>%
paste0('sigma', '[', . ,']')%>%sort
sigma <- matrix(sigma, neq, neq, byrow = TRUE)
sigma[lower.tri(sigma)] <- sort(sigma[upper.tri(sigma)])
#create Y
y <- paste0('eps[', 1:neq,']')
ypy <- paste0('[', y, ']', collapse = ',')
ypy <- paste0('Y : matrix(', ypy, ')')
#sigma
spy <- paste0(sapply(1:neq, function(i) paste0('[', paste0(sigma[i,], collapse = ', '), ']')),
collapse = ',')
spy <- paste0('Sigma : matrix(', spy, ')')
detS <- 'dets : ratsimp(determinant(Sigma))'
invS <- 'invs : ratsimp(invert(Sigma))'
expt <- 'expt : ratsimp(transpose(Y) . invs . Y)'
obj <- c(ypy, spy, detS, invS, expt)
#grind function in Maxima returns an object that can be mathematically evaluated
out <- c('print(new)', 'grind(dets)',  'print(new)', 'grind(invs)', 'print(new)', 'grind(expt)')
## Not run: 
rez <- wxMaxima(obj, out)

## End(Not run)