In the square-law equations you use (VDS - VDS,sat) while other books use

only VDS. Would you comment on the difference?

Yes, in Eq. (6.42) we derived, for saturation operation,

ID = [(KP/2)•(W/L)•(VGS - VTHN)2]•(1 + λ•(VDS - VDS,sat))

and for the triode region, Eq. (6.36)

ID = (KP)•(W/L)•[(VGS - VTHN)•VDS - V 2DS/2]

It's important to note that at the border between the triode and saturation

regions, where VDS = VDS,sat = VGS - VTHN, these equations are equal and

given by

ID = ID,sat = [(KP/2)•(W/L)•(VGS - VTHN)2]

Okay, now as you have indicated others remove the VDS,sat  and use simply

ID = [(KP/2)•(W/L)•(VGS - VTHN)2]•(1 + λVDS)

which is the implementation also used in the Level 1 SPICE model.

So here's the rub, if you don't use VDS,sat then there is a discontinuity

at the border of the triode and saturation regions when VDS = VDS,sat. To

avoid this discontinuity you have to also multiply the equation for

operation in the triode region by (1 + λVDS) or

ID = (KPW/L)•[(VGS - VTHN)• VDS - V 2DS/2]•(1 + λVDS)

This is often not mentioned or indicated by other authors, but, again it is what

is used in the level 1 MOSFET SPICE model.

See, also, the bottom of page 146 in the third edition.