The other approach is transforming to new coordinates similar to

http://www.mapleprimes.com/questions/100922-Pdsolve-Numeric-For-Multilayer-Heat-Conduction

transform-diffusion3.mw

yet here the problem lies with choosing the right transformation and i'm a little confused and unsure whether this approach is legit.

Also, how do i 3dplot my back-transformed results?

Still i didnt find the solution. I tried to find a way around maples reluctance to take contraints in the middle of the evaluated range (e.g. z(2000,t)=zeta(2000,t))

so i tried this

pde := [diff(z(x, t), t) = gDiffusion*10^6*(diff(z(x, t), x, x)), diff(zeta(x, t), t) = KDiffusion*10^6*(diff(zeta(x, t), x, x))];

ibc1 := zeta(x, 0) = .8, z(x, 0) = .8; #initial conditions
ibc2:=(D[1](zeta)(2000,t)=0, (D[1](z))(2000, t) = signum(cos((1/500)*t))/zF,  zeta(2000,t)=0.4, zeta(0,t)=z(0,t)) ; #constraints

sol := pdsolve(pde, {ibc1, ibc2}, numeric, time = t, range = 0 .. 2000, spacestep = 1);
sol:-plot3d(z(x, t), t = 0 .. 120)

like this zeta and z are only coupled at x=0 which seems to be an option

x for zeta and x for z need to be interpreted differently as it is the same coordinate in the calculation whereas in reality it is not...

still the problem is ill-defined and i am not happy with the "solution"