Creating a new varistor using SILVACO TCAD

hi.. i am new to Silvaco. I am facing difficulty in plotting the final graph. I am trying to create a simple structure of a varistor which is almost similar to diode where the base is TiO2 and dopants used are SiO2, Bi2O3, Ta2O5 and WO3 which will be varied for different experiments. The size of the varistor is 5mm.
I would also like to know how can i remove Titanium from the structure.
Besides that, I would really like to know how to create an IV graph for this structure. Below is what I have done so far. (im not sure if what I have done is correct)
thanking u in advance for your help....

go athena
#set mesh
#deposit TiO2

line x loc=0.00 spac=0.2
line x loc=1 spac=0.1
line x loc=1.1 spac=0.02
line x loc=2 spac=0.25


#
line y loc=0.0 spac=0.02
line y loc=0.2 spac=0.1
line y loc=0.4 spac=0.02
line y loc=2 spac=0.5


initial titanium
orientation=100 c.phos=1e14 space.mul=2

#pwell formationincluding masking off of the nwell
difuss time=30 temp=1000 dryo2 press=1.00 hcl=3

# Well oxidation
diffus time = 20 temp =970 dryo2
extract name = "Well Oxidation 200A "thickness material= "TiO~2"
mat.occno=1 x.val=0.05385
struct outfile= set2p-well_oxide.str

deposit material =TiO2 thick=0.5 divisions=2
struct outfile=set2mesh2.str


deposit material= SiO2 thick=0.02 divisions=3
struct outfile=set2mesh3.str


deposit material= WO3 thick=0.02
struct outfile= set2mesh4.str

deposit material = Ta2O5 thick=0.02
struct outfile= set2mesh5.str

deposit material= Bi2O3 thick=0.02
struct outfile=samplez.str


etch poly left p1.x=10
struct outfile=set2mesh6.str
electrode name= anode x=0.0
electrode name=cathode backside


struct outfile=set2meshtrial3.str

tonyplot -st set2meshtrial3.str
#quit

go atlas

mesh infile=set2meshtrial3.str



#region number=1 material=Bi2O3 y.max=-0.58 y.min=-0.56
#region number=2 material=Ta2O5 y.max=-0.56 y.min=-0.54
#region number=3 material=WO3 y.max=-0.54 y.min=-0.52
#region number=4 material=SiO2 y.max=-0.52 y.min=-0.5
#region number=5 material=TiO2 y.max=-0.5 y.min=0
#region number=6 material=Titanium y.max=0

material material = Bi2O3 user.group=semiconductor user.default=TiO2

#SRH
material material=Bi2O3 permittivity=30
material material=Bi2O3 eg300=2.7
material material=Bi2O3 mun=1.6 mup=1.6



doping uniform concentration=1e16 n.type region=2
doping uniform concentration=1e16 n.type region=3
doping uniform concentration=1e16 n.type region=4
doping uniform concentration=1e16 n.type region=5
doping uniform concentration=1e16 p.type region=6


#material material = Bi2O3 user.group=semiconductor user.default=TiO2
#dosextract eg300=2.9

#material material=Bi2O3 Young.m=3.2e10 Poiss.r=0.33
#material material=Bi2O3 components=2 density=8.9 abund.1=1 abund.2=2
#material material= Bi2O3 tcon.const=0.08
#material material=Bi2O3 visc.0= 5.154e8

method newton

# Bias the drain
solve vdrain=0.1

# Ramp the gate
log outf=trial2_1.log master
#solve vgate=0 vstep=0.25 vfinal=3.0 name=gate
#save outf=trial2_1.str


#contact name=gate n.TiO2
#interface qf=3e10


# extract device parameters
extract name="nvt" (xintercept(maxslope(curve(abs(v."gate"),abs(i."drain")))) \
- abs(ave(v."drain"))/2.0)
extract name="nbeta" slope(maxslope(curve(abs(v."gate"),abs(i."drain")))) \
* (1.0/abs(ave(v."drain")))
extract name="ntheta" ((max(abs(v."drain")) * $"nbeta")/max(abs(i."drain"))) \
- (1.0 / (max(abs(v."gate")) - ($"nvt")))


#Ramp the gate
log outf=trial3_1.log master
solve vgate=0 vstep=0.25 vfinal=3.0 name=gate

save outf=trial3.str

# plot results
tonyplot trial3_1.log -set trial3_1_log.set


quit