qs : [rho,m,E];

depends(qs,[t,x]);

"derivatives";

diff_t_rhs : [-m, (1-%gamma)*(E-1/2*m^2/rho)-m^2/rho, 1/2*(%gamma-1)*m^3/rho^2-%gamma*E*m/rho];

"evolution equations";

eqn_diff_t : makelist(diff(qs[i],t)='diff(diff_t_rhs[i],x),i,length(qs));

"discrete derivative";

"alpha is the +1 rhs coefficient, beta is the -1 coefficient, feel free to expand/exchange these for higher order spatial derivatives";

disc_subst_r : [rho=rho_r,m=m_r,E=E_r];

disc_subst_l : [rho=rho_l,m=m_l,E=E_l];

eqn_diff_t_disc : map(lambda([expr],alpha*subst(disc_subst_r,expr)+beta*subst(disc_subst_l,expr)),diff_t_rhs);

neighborVars : [rho_r,m_r,E_r,rho_l,m_l,E_l];

"first derivatives";

eqn_diff_t_q_disc : ratsimp(apply(matrix,map(lambda([expr], map(lambda([q],diff(expr,q)),neighborVars) ), eqn_diff_t_disc)));

"second derivative diff_t rho";

ratsimp(transpose(makelist(
apply(matrix,map( lambda([q], makelist(diff(eqn_diff_t_q_disc[i,j],q),j,matrix_size(eqn_diff_t_q_disc)[2])),neighborVars)),
i,matrix_size(eqn_diff_t_q_disc)[1])));