The main program sets up geometry of solution domain, and boundary conditions and initial conditions.

• The solution domain can be either 1-D or 2-D.
• Boundary conditions are applied on half-volumes, such that nodes are on boundaries, as described in Patankar.

Second m-file sets up and solves equations in a time-stepping loop. Volume sizes can be spatially variable, but I have not tested this feature - (i.e. "caveat emptor" :-)

• For illustration only, the code includes a parameter "theta" controlling the balance between past and future times
• theta = 1 fully implicit
• theta = 0 explicit scheme
• If using this code for a real problem,, clean out the theta code and uncomment the hard-wired fully implicit lines.
• visualization is still totlaly clunky

Solver solve.m (same old solver) solves the resulting system of linear equations at every time step.