All published worksheets from http://sagenb.org
# profit on a dominant assurance contract K,V = var("K,V") N = 100 sigma, mean = 1/12, 0.5
def G(x): return (1+erf((x-mean)/(sigma*sqrt(2.0))))/2
G(0.2)
def z(V,K): return V*K*binomial(N,floor(K))*(1-G(V))^K*G(V)^(N-K)
z(0.5,40)
parametric_plot3d([K, V, z], (K,0, 100), (V, 0, 1), plot_points=[101,51], mesh=True)