Subscript[u,1]

n=4
Δz=L/n;
eq1=Table[(cD/((Δz)^2))(Subscript[u,j+1][t]+Subscript[u,j-1][t]-2Subscript[u,j][t])==D[Subscript[u,j][t],t]+g Subscript[u,j][t],{j,1,n-1}]

cf={Subscript[u,0][t]->Subscript[u,1][t],Subscript[u,n][t]->Subscript[u,n-1][t]}
eq2=eq1/.cf

Plot[1/2(Tanh[((z-a)/s)]+1)/.{a->0.3,s->0.01},{z,0,1}]

Plot[(1/2)(Tanh[((z-a)/s)]-Tanh[((z-b)/s)])/.{a->0.3,b->0.7,s->0.001},{z,0,1}]

U0[z_]:=(u0/2)(Tanh[((z-a)/s)]-Tanh[((z-b)/s)])
ci=Table[Subscript[u,j][0]==U0[Δz j],{j,1,n-1}]

difeq=Join[eq2,ci]

ndifeq=difeq/.{u0->1,L->1,a->0.3,b->0.7,cD->0.1,g->0.5,s->0.005}

var=Table[Subscript[u,j][t],{j,1,n-1}]

(*No se olviden de matar el kernel aqui *)
n=100;
U0[z_]:=(u0/2)(Tanh[((z-a)/s)]-Tanh[((z-b)/s)])

Δz=L/n;
eq1=Table[(cD/((Δz)^2))(1+(0/2)(Tanh[((j Δz-0.5)/0.01)]+1))(Subscript[u,j+1][t]+Subscript[u,j-1][t]-2Subscript[u,j][t])==D[Subscript[u,j][t],t]-g,{j,1,n-1}];
cf={Subscript[u,0][t]->Subscript[u,1][t],Subscript[u,n][t]->Subscript[u,n-1][t]};
eq2=eq1/.cf;
ci=Table[Subscript[u,j][0]==U0[Δz j],{j,0,n}];
difeq=Join[eq2,ci];
conds={u0->1,L->1,a->0.3,b->0.7,cD->0.1,g->0.00,s->0.005};
ndifeq=difeq/.conds;
var=Table[Subscript[u,j][t],{j,1,n-1}];
sol=NDSolve[ndifeq,var,{t,0,200}];

listplot=Table[{(j Δz/.conds),(Subscript[u,j][t]/.cf/.sol[[1]]/.{t->1})},{j,0,n}];

ListPlot[listplot,Joined->True]

L

listplot=Table[{(j Δz/.conds),(Subscript[u,j][t]/.cf/.sol[[1]]/.{t->0.10})},{j,0,n}];
ListPlot[listplot,Joined->True,PlotRange->{-0.05,1.01}]

tf=1.0;
zm=1.0;
plot={};
ProgressIndicator[Dynamic[prog]]
Do[
prog=vz/zm;
Do[
list=Table[{(j Δz/.conds),(Subscript[u,j][t]/.cf/.sol[[1]]/.{t->vt})},{j,0,n}];
vu=Interpolation[list,vz];
AppendTo[plot,{vt,vz,vu}];
,{vt,0,tf,tf/40}];
,{vz,0,zm,zm/40}]

plot[[2]]

ListPlot3D[plot,PlotRange->All]

ListPlot3D[plot,PlotRange->All]

ListDensityPlot[plot,PlotRange->All,ColorFunction->"AvocadoColors"]