Implement the explicit numerical scheme for 1-D diffusion of a concentration u(x,t) 
in an interval [a, b], during time [0, tmax], with boundary conditions specified below.

Your code should accept as inputs: MM, tmax, dtout, factor, a, b, D ,
(then dx=1/MM, M=(b-a)*MM), and the time-step dt should be set as:
	dt = factor * dtEXPL 
The code should output the profile:   x(i)   U(i)  (i=0,...,M+1)   initially 
(at time=0), at desired times tout (every dtout), and also at time tmax, 
in a file for plotting, i.e. in 2 columns, like (see below):
        # Profile at time: ???? , nsteps= ????
          x(0)             U(0)
          ....		  ....
          x(M+1)	U(M+1)

An explicitly solvable problem on which you can debug your code
is the following diffusion problem in a semi-infinite interval:
        ut = D uxx , 0 < x < Inf , t > 0
        u(x,0) = 0,   u(0,t) = 1
whose exact solution is:   u(x,t) = ERFC( 0.5 x / SQRT(D t) ),
with ERFC() = 1 - ERF() the Complementary ERROR function. 
Verify (by hand) that this u(x,t) satisfies the problem.

To simulate this problem on a finite interval [0,b], we simply impose 
the exact solution at x=b, i.e. use the boundary conditions:
        U(0) = 1.0 ,   U(M+1) = ERFC( 0.5*b / SQRT(D*time) )
Thus, the exact solution and the numerical solution solve the same IBVP problem 
and we can compare the two to calculate the maximum error at any time tn:
        ERRn = MAX{ |U(i) - u(xi, tn)| : i=0,...,M+1 }.

Something like the following routine can be called at each time-step
to find the exact solution and the error: 

      subroutine COMPARE( M, D, x, time, U, uEXACT,  ERR )
!! Compare with exact solution: u(x,t) = ERFC(.5*x/SQRT(D*t)) 
!! of:     u_t = D u_xx , 0<x<Inf , t>0 ; u(x,0)=0 ; u(0,t)=1
      implicit real*8 (a-h, o-z)
      dimension x(0:M+1), U(0:M+1), uEXACT(0:M+1)
                ERR = 0.0
          DO i = 0, M+1
                   arg = 0.5 * x(i) / SQRT( D*time )
             uEXACT(i) = ERFC( arg )	!! = 1 - ERF(arg)
                  ERRi = ABS( U(i) - uEXACT(i) )
                  ERR  = MAX( ERRi, ERR )
          ENDDO
      END


1. Debug using:  a=0., b=1., D=0.1, MM=8, tmax=3., dtout=1., factor=0.98 
    With only 8 nodes (MM=8), the agreement won't be very good.
    Nevertheless, try  factor=0.90 to see how it improves.

2. Try MM=16 to see how it improves;  the plotted curves should be indistinguishable, 
    and the error at tmax=3 should be less than 1%.

3. Try MM=32 and MM=64 to see further improvement in the error.
    Make runs with factor = 0.90 , 0.99 to see the effect.
    Now you have a debugged diffusion code !

4. To see how easy it is to go unstable, make runs with factor=1.0,
    and then with factor=1.01 and 1.1. Look at the plot and the error.

Submit the ERRor  at tmax from each of 1., 2., 3., 4. , and your code.


                            Remarks - Hints
a. Only x, F, U, uEXACT  need to be indexed as (spatial) arrays. 
   The natural indexing for M control volumes is:
	a=x(0),  x(1), ..., x(M),  x(M+1)=b
	[i=1,...,M are the nodes, and i=0, i=M+1 are the boundary faces]
	  U(0),  U(1), ..., U(M),  U(M+1)
	F(1),  F(2), ..., F(M),  F(M+1)		( F(i)=Fi-1/2 )
	[i=2,...,M are the internal fluxes, and i=1, i=M+1 are the boundary fluxes]
   [For Matlab coding, indices must be shifted to start with 1, unfortunately]

b. It is highly advisable to use modular programming, with different
   tasks in separate subprograms (or at least organize your code), like:
	MESH	construct the mesh: x(i), i=0,...,M+1
	INIT	set initial values for U(i), i=0,...,M+1
	FLUX	compute fluxes F(i), i=1,...,M+1  ( F(i)=Fi-1/2 )
	PDE	update U(i), i=0,...,M+1
	COMPARE	compute exact solution and error
	OUTPUT	print output

   The structure of the main program will be something like this

	declare variables, dimension arrays 
	read in the data from dat file (see below)
	    call MESH( ... ) 	
	!............ set the timestep ...........
	  dtEXPL = dx*dx/(2*D)	! max timestep for stability
	  dt     = factor*dtEXPL	! factor<1 increases stability
	!............... initialize ...........
	    call INIT( ... )
	  nsteps = 0
	    time = 0.0
	    tout = MAX( dtout, dt )  
	!............... begin time-stepping ...............
        DO nsteps = 1, Nmax
	      time = time + dt	!!or: = nsteps*dt
	      call FLUX( ... )
	      call PDE ( ... )
	      call COMPARE( ... )  
	         
	    IF( time >= tout ) THEN
		call OUTPUT( ... )
		tout = tout + dtout 
	    ENDIF
	    IF( time >= tmax ) THEN
	       print*, ' DONE, at time= ', time, '  nsteps=', nsteps
	       print*, '       max error =', ERR
	 	STOP
	    ENDIF
        ENDDO !.......... end of time-stepping .............
     	print*,'... out of timesteps: need bigger Nmax'
	END

c. The input data should be read in at run time so they can be changed easily for each run. 
    They are too many to type each time, so the best way is to read them from a dat file.
    e.g. here is a convenient arrangement (in a file dat):
	8    3.   1.                  : MM, tmax, dtout
	0.98    0.1  0.0   1.0    : factor,  D,  a, b 
  [ In matlab, can use a script m-file  dat.m:
      M = 8 ;  tmax=3. ;  dtout=1. ;  Nmax=1000; factor=0.98 ;  D = 0.1 ;  a=0.0 ; b=1.0  
     which the code reads with:  dat  ]

   These are good values to start with. On unix, you execute with:
	a.out < dat        (output to screen), or
	a.out < dat > out  (output to file out)

d. The OUTPUT routine, called at each tout (every dtout)
   after COMPARE(...), should record the time and nsteps
   of the profile, so it should be something like:

      print"( '# Profile at time:', f9.4,'  nsteps=', I6 )",time,nsteps
      print"( '# Error up to this time:', e15.6 )", ERR
      DO i = 0, M+1
           print"( f12.4, g16.8, g16.8 )", x(i), U(i), uEXACT(i)
      ENDDO
      print"()"       !! blank line for gnuplot to break plot

   Send your output to a file (a.out < dat > out) and can use gnuplot 
   to see the numerical and exact solutions:
       gnuplot> plot "out" u 1:2 w lines lt 1 lw 2, "out" u 1:3 w lines lt 3 lw 2

   To get the curves only initially and at tmax=3.0, run with dtout=3.

e. For Fortran programmers: in double precision, you must use DERFC 
    for the (double precision) complementary error function.