;lecture_diffusion_task1c
;linear part of the solution
;Does not depend on time !!
linpart=temp0+((temp1-temp0)/Lx)*x
fxhelp=fx0-1.0*linpart
;
;Compute Fourier components for -linpart
maxk=nx-1; how may coeficients?
ak_vec=fltarr(maxk+1) ; fourier coefficients
for k=1,maxk do begin
help1=total(fxhelp*sin(k*!pi*x/Lx))
ak_vec[k]=(2.0*dx/Lx)*help1
endfor
;
;Compute and display time dependent solution
;
erase
;nt=1 ; if nt=1 then we just show rho(x,0)
for it=0,nt-1 do begin; do begin; time evolution
t=tvec[it]
rho=linpart
;Compute solution for time t
for k=1,n_elements(ak_vec)-1 do begin
ak=ak_vec[k]
dummy=ak*sin(k*!pi*x/Lx)
dummy2=exp(-(k^2*!pi^2*D*t)/Lx^2)
rho=rho+dummy*dummy2
endfor
plot,x,rho,yrange=1.0*[-0.2,1.2*A],linestyle=0,thick=2,$
xtitle='x',ytitle='rho',title='1D Diffusion equation',$
charsize=1.5
oplot,x,fx0,linestyle=1
xyouts,2,-0.3,'t='+string(strcompress(t)),size=2
xyouts,0,temp0-0.1,'O',size=3,charthick=2,alignment=0.5
xyouts,Lx,temp1-0.1,'O',size=3,charthick=2,alignment=0.5
wait,0.2
endfor