My code is here:
* MEMRISTOR
* Ron, Roff - Resistance in ON / OFF States
* Rinit - Resistance at T=0
* D - Width of the thin film
* uv - Migration coefficient
* p - Parameter of the WINDOW-function
* for modeling nonlinear boundary conditions
* x - W/D Ratio, W is the actual width
* of the doped area (from 0 to D)
*
.SUBCKT memristor Plus Minus
+ Ron=100 Roff=16K Rinit=11K D=10N uv=10F p=10
* DIFFERENTIAL EQUATION MODELING *
.PARAM f(x,p)= '1-pow((2*x-1),(2*p))'
Gx 0 x value='I(Emem)*uv*Ron/pow(D,2)*f(V(x),p)'
Cx x 0 1 IC='(Roff-Rinit)/(Roff-Ron)'
Raux x 0 1T
* RESISTIVE PORT OF THE MEMRISTOR *
Emem plus aux value='-I(Emem)*V(x)*(Roff-Ron)'
Roff aux minus Roff
*Flux computation*
*Eflux flux 0 value={SDT(V(plus,minus))}
*Charge computation*
*Echarge charge 0 value={SDT(I(Emem))}
* WINDOW FUNCTIONS
* FOR NONLINEAR DRIFT MODELING *
*window function, according to Joglekar
*.func f(x,p)={1-(2*x-1)^(2*p)}
*proposed window function
*;.func f(x,i,p)={1-(x-stp(-i))^(2*p)}
.ENDS memristor
Xmemrist aa 0 memristor Ron=100 Roff=16K Rinit=11K D=10N uv=10F p=10
Vtest aa 0 SIN(0 1.2V 1 0 0 0)
.tran 1m 1s
.plot tran V(aa) I(Vtest)
.option list node post=2 ingold=2 runlvl=0
.end
and my results are like this: