(*  Begin with the assumption that velocity is constant                  *)
(*                                                                       *)
v[x,t]=v0
(*                                                                       *)
(* 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).  To reproduce   *)
(*  something that looks like a differential equation we divide dm by dt.*)
(*                                                                       *)
diffdm= Expand[dm/dt]
(*                                                                       *)
(*  Our next goal is to manipulate this mess, so that we can get some-   *)
(*  thing that looks like the diffusion term in an advection diffusion   *)
(*  equation.                                                            *)
(*                                                                       *)
(*  Mathematica isn't to bright about rules related to derivatives in    *)
(*  general or the relations forced by mass concervation in particular   *)
(*  so we must force some substitutions.  In doing so the starting point *)
(*  is the 1-D differential equation for conservation of mass.  The      *)
(*  first substitution comes from taking the partial derivative with     *)
(*  respect to time of the entire mass equation, giving a way to         *)
(*  eliminate the second derivatives with respect to time.               *)
(*                                                                       *)
temp1= diffdm /. D[rho[x,t],{t,2}]->D[-rho[x,t]*v[x,t],x,t]
(*                                                                       *)
(*  The above line simple says "Replace (/.) all occurances of the 2nd   *)
(*  derivative with respect to time of density with the partial with     *)
(*  respect to time of minus the mass flux term.  The replacement is     *)
(*  done on the contents of diffdm and placed in temp1.                  *)
(*                                                                       *)
(*  Next we get rid of all cross derivatives by looking a the derivative *)
(*  with respect to x of the mass equation.                              *)
(*                                                                       *)
temp2 = temp1 /. D[rho[x,t],x,t]->D[-rho[x,t]*v[x,t],{x,2}]
(*                                                                       *)
(*   Finally we substitute in the time derivative of density.  This      *)
(*   cleans up a couple odds and ends and leaves only terms that the     *)
(*   difference scheme effectly adds to the normal mass conservation     *)
(*   equation.                                                           *)
(*                                                                       *)
temp3 = temp2  //. D[rho[x,t],t]->D[-rho[x,t]*v[x,t],x]
(*                                                                       *)
(*     "//." says to keep making the substitution until nothing changes  *)
(*  We reorganize temp3 and put it in the Mathematica variable "resid"   *)
(*  to represent the numerical residual to the correct mass equation.    *)
(*                                                                       *)
resid = Factor[temp3]