If you haven't done so already, modify your program from HW2 to use 2nd order accurate surface flux terms (fit a quadratic curve to the surface temperature and the two volume centers in from the surface).  Note that I'm requiring a second order evaluation of the temperature derivative at the surface in your flux term.  I am not requiring that you get a conduction equation in the surface volumes that is second order accurate.   If your HW2 implementation was first order, provide me with a revised set of equations for the 9 volume types. 

Run the HW2 problem with 3x3, 9x9, and 27x27 volumes.  Apply Richardson Extrapolation analysis to temperatures at the center point in the rod, and points located midway along one x-facing surface and midway along one y-facing surface.  Note that the last two evaluations are surface not volume centered temperatures and you will need to use an equation consistent with your second order surface flux evaluation to get them.  At each of these three points in space give me the predicted order of accuracy, the predicted error on the finest (27x27) mesh, and predicted error on the 9x9 mesh.

You will probably need more decimal digits of precision for your temperatures than currently available in the file steady.out.   Just after the line:

21  FORMAT(10F8.2)

in subroutine edit, add the following lines:

    WRITE(11,*)
    WRITE(11,*) ' Middle of y facing edge T = ', tys(nx/2+1,1)
    WRITE(11,*) ' Middle of x facing edge T = ', txs(ny/2+1,1)
    WRITE(11,*) ' Center Point in the rod T = ', t(nx/2+1,ny/2+1)