c*************************************************************
c*
c*  C-N and Implicit methods for solution of heat conduction problem...
c*
c*   diff(u,z,z) = diff(u,t) + f(x,t)  for 0<x<xR  and 0<t<tmax
c*
c*   where u(x,0)=g(x), u(0,t)=uL, and u(L,t)=uR
c*
c*************************************************************
      implicit real*8(a-h,o-z)
      real*8 lam, L
      parameter(L=1.0,uL=0.0,uR=0.0,nx=100,nt=1,tmax=0.001)
      parameter(a=0.25,b=0.8, pi=3.141592654)
      dimension x(0:nx+1),t(0:nt),ucn(0:nx+1),uip(0:nx+1)
      dimension w(nx),v(nx),z(nx),uex(0:nx+1)
      f(x,t)=0.0
c      g(x)=sin(2.0*pi*x)
c      g(x)=1
      g(x)=step(x-a)-step(x-b)
      icase=2
      dx=L/(nx+1)
      dt=tmax/nt
      lam=dt/(dx*dx)	
      t(0)=0.0
      open(unit=6,file="data")
      write(6,106)
  106 format('x-axis,',' Exact Solution,',' Implicit Method,',
     &' Crank-Nicolson Method')
c  
c   Crank-Nicolson
c  
      do 2 ix=0,nx+1
        x(ix)=dx*ix
    2   ucn(ix)=g(x(ix))
      xa=-2.0*(1.0+lam)
      do 20 it=1,nt
        t(it)=dt*it
	w(1) = xa
        v(1) = lam/w(1)
	ff=dt*(f(x(1),t(it))+f(x(1),t(it-1)))
	uu=-lam*ucn(2)-2.0*(1.0-lam)*ucn(1)-lam*ucn(0)
        rr=uu+ff
	z(1)=rr/w(1)
        do 12  jj = 2, nx
	   w(jj)=xa-lam*v(jj-1)
	   v(jj) = lam/w(jj)
	   ff=dt*(f(x(jj),t(it))+f(x(jj),t(it-1)))
	   uu=-lam*ucn(jj+1)-2.0*(1.0-lam)*ucn(jj)-lam*ucn(jj-1)
           rr=uu+ff
   12      z(jj) = (rr-lam*z(jj-1))/w(jj)
	ucn(nx) = z(nx)
	do 16  jj = nx - 1, 1, -1
   16      ucn(jj) = z(jj)-v(jj)*ucn(jj+1)
        ucn(0)=uL
        ucn(nx+1)=uR
   20   continue
c  
c   Implicit
c  
      xa=-1.0-2.0*lam
      do 22 ix=0,nx+1
   22   uip(ix)=g(x(ix))
      do 30 it=1,nt
	w(1) = xa
        v(1) = lam/w(1)
        rr=-uip(1)+dt*f(x(1),t(it))
	z(1)=rr/w(1)
        do 24  jj = 2, nx
	   w(jj)=xa-lam*v(jj-1)
	   v(jj) = lam/w(jj)
           rr=-uip(jj)+dt*f(x(jj),t(ij))
   24      z(jj) = (rr-lam*z(jj-1))/w(jj)
	uip(nx) = z(nx)
	do 26  jj = nx - 1, 1, -1
   26      uip(jj) = z(jj) -v(jj)*uip(jj+1)
        uip(0)=uL
        uip(nx+1)=uR
   30   continue
      if(icase.eq.1)  call exact1(uex,x,tmax,nx,pi)
      if(icase.eq.2)  call exact2(uex,x,tmax,nx,pi,a,b)
      open(unit=6,file='data')
      write(6,100)  (x(ix),uex(ix),uip(ix),ucn(ix),ix=0,nx+1)
  100 format(e10.4,',',e12.4,',',e12.4,',',e12.4)
      close(unit=6)
      end

      function step(x)
      implicit real*8(a-h,o-z)
      step=0.0
      if(x.ge.0.0) step=1.0
      return
      end

      subroutine exact1(uex,x,t,nx,pi)
      implicit real*8(a-h,o-z)
      dimension uex(0:nx+1),x(0:nx+1)
      do 10 ix=0,nx+1
   10   uex(ix)=exp(-4.0*pi*pi*t)*sin(2*pi*x(ix))
      return
      end
      
      subroutine exact2(uex,x,t,nx,pi,a,b)
      implicit real*8(a-h,o-z)
      dimension uex(0:nx+1),x(0:nx+1)
      do 10 ix=0,nx+1
        sum=0.0
        do 8 i=1,100
          c=2.0*(cos(a*i*pi)-cos(b*i*pi))/(i*pi)
    8     sum=sum+c*exp(-i*i*pi*pi*t)*sin(i*pi*x(ix))
   10   uex(ix)=sum
      return
      end