pro dftcs

; FTCS for the Diffusion Equation
; "Temperature" diffuses in time from an initial delta funtion spike

	dt = 2.0d-4		;time step
	N  = 41			;grid points
	L  = 1.d0		;system size (from -L/2 to L/2)
	kappa  = 1.d0		;diffusion coefficient

	h  = L/(N-1)		;step size
	coeff = kappa*dt/h^2

	print,' coeff = ',coeff

	tt = fltarr(N)*0.d0
	tt_new = tt
	tt(N/2) = 1/h		;delta function


	iplot = 0
	nstep = 250
	nplots = 25
	tt_plot = fltarr(N,nplots)*0.d0  ; holds 'temperature profile'
	xplot = findgen(N)*h - L/2       ; xvals for plot
	tplot = fltarr(nplots)           ; will hold tvals for plot


	plot_step = nstep/nplots

	for istep=1,nstep do begin
		tt_new(1:N-2) = tt(1:N-2) + $
			coeff*( tt(2:N-1) + tt(0:N-3) - 2*tt(1:N-2))
		tt = tt_new
		if ((istep mod plot_step) eq 0) then begin
			tt_plot(*,iplot) = tt
		        tplot(iplot) = istep*dt
			iplot = iplot + 1
			print,' No of steps completed = ',nstep
		endif
	endfor

	surface,tt_plot,xplot,tplot,xtitle='Distance',ytitle='Time',$
	        title='Diffusion of delta function spike',charsize=2
	return
	end