restart;
C0 := 10:
K  := 3.0e6:
d  := 5.0E-14:
hm := 5.0E-4:
L  := 0.001:

DE1 := diff(C(x, t), t) = d*diff(C(x,t), x, x);
IC := C(x, 0) = piecewise( x < 0, 0,
                           0 <=x and x <= L, C0,
                           L < x, 0
                         );
BC1 := D[1](C)(0, t) = 0;
BC2 := -d*D[1](C)(L, t) = hm * C(L, t) / K;
pds:=pdsolve( DE1, {IC, BC1, BC2},
              type=numeric, 
              spacestep=L/10000,
              timestep=10000);

T := 60*60*24*100;

p1:=pds:-plot3d(t=0..T, x=0..L);
plots:-display(p1);
p2:=pds:-plot(x=L, t=0..T, color="Brown");
#plots:-display(p2);

v1:=pds:-value(x=L, t=0..T,output=listprocedure);
CatL:=rhs(v1[3]);
v1Vals:=[ seq( [j, CatL(j)], j=0..T,T/10)];

CFit1 := CurveFitting:-Spline([seq([j, CatL(j)], j = 0 .. T, (1/1000)*T)], t):
CFit:=unapply(CFit1, t):

p3:= plot(CFit, 0..T, color="Blue"):
plots:-display( [p2,p3]);

seq([j, fdiff(CFit(x),x,x=j)], j=0..L,L/10);

#3. Check flux;
FluxVals1:=[seq([j, -d*fdiff(CFit(t),t, t=j)], j=0..T,T/10)];
p6:=plot(-d*diff(CFit(t),t), t=0..T, color="Brown");
#p6b:=plot(-d*fdiff(CFit(t),t, t=tau), tau=0..T, color="Red");
plots:-display(p6);

Flux:=proc (t) options operator, arrow; hm/K*CatL(t) end proc;
FluxVals2:=[seq([j, Flux(j)], j=0..T,T/10)];
p7:=plot(Flux(t),t=0..T, color="Purple");
plots:-display(p7);

#Mass emitted;
CatT:=pds:-value(t=T, output=listprocedure);
func_CatT:=rhs(CatT[3]);
#plot(func_CatT, 0..L):
IntMdot := C0*L-int(func_CatT(x), x=0..L);
AverageFlux := IntMdot/T;