U_perturbed_def : U = U_o + %epsilon * U_hat * exp((y*m + z*k)*%i);

/* TODO rewrite the above eqn per-variable but with all prims on the lhs */"";

E_from_P : E = %rho * (ux^2 + uy^2 + uz^2)/2 + P / (%gamma - 1);

U = transpose(matrix([%rho, %rho * ux, %rho * uy, %rho * uz, E]))$
subst([E_from_P],%)$
U_def : %;

Uo_def : subst([%rho=rho_o, ux=ux_o, uy=uy_o, uz=uz_o, P=P_o, U=U_o], U_def);

Uhat_def : subst([%rho=rho_hat, ux=ux_hat, uy=uy_hat, uz=uz_hat, P=P_hat, U=U_hat], U_def);

W_def : W = transpose(matrix([%rho, ux, uy, uz, P]));

Wo_def : W_o = transpose(matrix([rho_o, ux_o, uy_o, uz_o, P_o]));

What_def : W_hat = transpose(matrix([rho_hat, ux_hat, uy_hat, uz_hat, P_hat]));

jacobian(transpose(U)[1],transpose(W)[1])$
subst([U_def, W_def], %)$
ev(%,transpose,jacobian)$
dU_dW_def : %;

invert(dU_dW_def)$
ratsimp(%)$
dW_dU_def : %;

us : [ux,uy,uz];

makelist( F[i] = transpose(matrix([
    %rho * us[i],
    %rho * us[i] * ux + if i=1 then P else 0,
    %rho * us[i] * uy + if i=2 then P else 0,
    %rho * us[i] * uz + if i=3 then P else 0,
    us[i] * (E + P)
])), i, 1, 3)$
subst([E_from_P],%)$
ratsimp(%)$
Fdefs : %;

xs : [x,y,z];

eqn1 : 'diff(U,t) + sum('diff(F[i],xs[i]), i, 1, 3) = 0;

dF_dx_expanded_def : makelist(
    'diff(F[i],xs[i]) = 'diff(F[i],W) . 'diff(W,U) . 'diff(U,xs[i])
,i,1,3);

eqn2 : subst(dF_dx_expanded_def, eqn1);

makelist(F[i],i,1,3)$
subst(Fdefs,%)$
makelist( 'diff(F[i],W) = jacobian( transpose(%[i])[1], transpose(subst([W_def],W))[1] ), i,1,3)$
dF_dW_defs : %;

makelist( 'diff(F[i],U) = ratsimp(rhs(dF_dW_defs[i]) . dW_dU_def) , i,1,3);

/* maxima sucks, because subst *automatically simplifies* after substituting, and sublis doesn't ... but *won't work* on subscript variables, nor do depends() work on subscript variables */

U_perturbed_def$
sublis([Uhat_def, Uo_def, U_def],%)$
makelist(lhs(%)[i,1]=rhs(%)[i,1],i,1,length(lhs(%)))$
U_perturbed_defs : %;

eqn2$
subst(['diff(W,U) = dW_dU_def], %)$
subst(dF_dW_defs, %)$
subst(U_def, %)$
subst(U_perturbed_defs, %)$
eqn3 : %;

append(transpose(rhs(Wo_def))[1], transpose(rhs(What_def))[1])$
depends(%,[t,x,y,z]);

eqn3$
ev(%,diff)$
makelist(lhs(%)[i,1]=0,i,1,length(lhs(%)));