linear(vds,vgsvth) := 2*vgsvth*vds-vds^2;
saturation(vds,vgsvth) := vgsvth^2;
draincurrent(vds,vgsvth) := if (vds>vgsvth) then saturation(vds,vgsvth) else linear(vds,vgsvth);
limit(vds) := vds^2;
load(draw);
draw2d(terminal=svg,file_name="IvsV_mosfet",dimensions=[480,320],
  color=black,font="Times",nticks=75,proportional_axes=none,font_size=12,
  user_preamble="set termoption enhanced;set mxtics; set mytics; set key bottom right",
  xtics_axis=false,ytics_axis=false,grid=true,points_joined=true,yrange=[0,50],
  xaxis=false,yaxis=true,axis_top=true,axis_bottom=true,axis_left=true,axis_right=true,
  line_width=2,xlabel="Drain to source voltage [V]",
  ylabel="Drain current [arbitary unit]",
  label(["linear region",2,40], ["saturation region",6.5,20]),
  label_alignment=right,label(["V_{GS}-V_{TH}=7 V",9.8,47],["6 V",9.8,37.5],
    ["5 V",9.8,26.5],["4 V",9.8,17.4],
    ["3 V",9.8,10.5],["2 V",9.8,5.5],["1 V",9.8,2.5]),
  color=blue,key="",explicit(draincurrent(x,1),x,0,10),explicit(draincurrent(x,2),x,0,10),
  explicit(draincurrent(x,3),x,0,10),explicit(draincurrent(x,4),x,0,10),
  explicit(draincurrent(x,5),x,0,10),explicit(draincurrent(x,6),x,0,10),
  explicit(draincurrent(x,7),x,0,10),
  color=red,line_width=1,line_type=dots,nticks=10,explicit(limit(x),x,0,10)
 );