library(gstat)
library(sp)

p = SpatialPoints(cbind(c(0,0,3),c(0,2,1)))
p$z = c(10,7,3)
krige(z~1, p, SpatialPoints(cbind(1,1)), vgm(psill=3, nugget=1, range=6, "Exp"))

#Results
[using ordinary kriging]
  coordinates  var1.pred var1.var
1      (1, 1)  6.813431  2.002507

σ²ε = sill - [w1 ... wn λ] [C10]
                           |...|
                           |Cn0|
                           [1  ]

C(h) = σ² = 4, if |h| = 0
C(h) = (σ² - a)exp(-3|h|/r) = (4 - 1)exp(-3|h|/6), if |h| > 0

C = function(h){3*exp(-3*h/6)}

VC = function(h){4-C(h)}

VCgstat = vgm(psill=3, nugget=1, range=6, "Exp")
plot(variogramLine(VCgstat,maxdist=50), type="l")
lines(0:50, VC(0:50), col="red")

VCgstat2 = vgm(psill=3, nugget=1, range=2, "Exp")

> krige(z~1, p, SpatialPoints(cbind(1,1)), vgm(psill=3, nugget=1, range=2, "Exp"))
[using ordinary kriging]
  coordinates var1.pred var1.var
1      (1, 1)  6.887426 3.136307