Math578 - Alexiades
                    LAB 3
        1D Diffusion - Comparison with exact solution
 
Consider a (1-dimensional) diffusion process, with (constant) 
diffusivity D > 0, for the concentration u(x,t) in the 
semi-infinite interval 0 ≤ x < ∞, during some time  0 ≤ t ≤ tend.  
Initialy, the concentration is uniform: uinit(x)=0, 
and the endpoint at x=0 is raised and held at concentration u0(t)=1. 
This problem admits the explicit similarity solution:

		u(x,t) = 1 − ERF( x / (2 SQRT(D t)) ) ,

with ERF(z) the error function.  
To test your diffusion code against this exact solution on an interval [0,b], 
need to impose the exact solution at x=b:

	U(M+1) = 1.0−ERF( b / (2*SQRT(D*t)) ) == ERFC( b / (2*SQRT(D*t)) ) .

So here both BCs are Dirichlet type, and the one at x=b is time-dependent...
think when time needs to be updated...

Take b=1.0, D=0.1, factor=0.96, and compare the numerical and exact solutions 
at tend=1, on a mesh with Δx=1/16 (MM=16) and a finer mesh with Δx=1/64 (MM=64).

The comparison should be done in a subprogram COMPARISON, with
minimal changes to the code (so that it can be turned off easily, by commenting out,  
on problems without exact solution)

It should print out the numerical solution U(i) at tend, the exact solution u_exact=u(xi, tend), 
and the error (their difference), i.e. print:
		xi     Ui      uEXACT     error
Also compute the worst absolute error ( maximum |error| over all i ).

You can plot both U and u_exact on the same plot in gnuplot:
    gnuplot>:  plot "out" u 1:2  lt 3 lw 1 , "out" u 1:3  lt 2 lw 2
(u 1:3  means using columns 1 and 3 as x and y respectively).

Prepare a text file "Lab3_code.txt" containing:
	1. max error for MM=16 
	2. max error for MM=64
	======================
	conclusions 
	======================
	your code
and a PDF file "Lab3_plots.pdf" containing the (annotated) plots.
If you can, combine them into a single PDF file "Lab3.pdf".
Submit on Canvas.
 
CAUTION for fortran programmers: The double precision error function in most fortran
libraries is "DERF", and ERFC is DERFC. It should be coded as DERFC( 1.d0 / (2*SQRT(D*t)) )