(* First we define a varable dm, that you should recognize as the left *)
(* hand side of a fully implicit mass conservation equation using a    *)
(* donor cell average and assuming a positive velocity                 *)
(*                                                                     *)
dm=
    rho[i,n+1]-rho[i,n]+dt/dx*(rho[i,n+1]*v[i+1/2,n+1]-rho[i-1,n+1]*
    v[i-1/2,n+1])
(*                                                                       *)
(*  Next write  the density at spatial cell center j and time level n as *)
(*  a Taylor expansion from cell center j at time level n+1              *)
(*  To permit proper use of derivatives in mathematica we will replace   *)
(*  the cell index j with the coordinate x and the time level index n+1  *)
(*  with t                                                               *)
(*                                                                       *)
(*  The Mathematica function D performs derivatives.  D[f[x,t],t] is the *)
(*  partial derivative of the function f with respect to t.  D[f[x,t],t,x] *)
(*  is the partial derivative with respect to t of the partial derivative  *)
(*  with respect to x of the function f.  D[f[x,t],{t,2}] is the second  *)
(*  derivative of f with respect to t.                                   *)
(*                                                                       *)
rho[i,n]  = rho[x,t]-dt*D[rho[x,t],t]+1/2*dt^2*D[ rho[x,t],{t,2}]
(*                                                                       *)
(*  Write  the density at spatial cell center j-1 and time level n+1 as  *)
(*  a Taylor expansion from cell center j at time level n+1              *)
(*                                                                       *)
rho[i-1,n+1] = rho[x,t]-dx*D[rho[x,t],x]+1/2*dx*dx*D[ rho[x,t],{x,2}]
(*                                                                       *)
(*  Write  the velocity at spatial cell face j+1/2 and time level n+1 as *)
(*  a Taylor expansion from cell center j at time level n+1              *)
(*                                                                       *)
v[i+1/2,n+1] = v[x,t]+1/2*dx*D[v[x,t],x]+1/8*dx*dx*D[ v[x,t],{x,2}]
(*                                                                       *)
(*  Write  the velocity at spatial cell face j-1/2 and time level n+1 as *)
(*  a Taylor expansion from cell center j at time level n+1              *)
(*                                                                       *)
v[i-1/2,n+1] = v[x,t]-1/2*dx*D[v[x,t],x]+1/8*dx*dx*D[ v[x,t],{x,2}]
(*                                                                       *)
(*  To keep notation straight we must redefine the density at cell j and *)
(*  time level n+1 in terms of x and t                                   *)
(*                                                                       *)
rho[i,n+1] = rho[x,t]
(*                                                                       *)
(*  Mathematica automatically substitutes all of the above replacement   *)
(*  equations into the expression for dm.  Now we tell Mathematica to    *)
(*  spread out all the terms in the expression (Expand) and regroup them *)
(*  as coefficients of powers of the cell length dx.  To reproduce       *)
(*  something that looks like a differential equation we divide dm by dt.*)
(*                                                                       *)
diffdm= Collect[Expand[dm/dt],dx]