function u=MLP (n,m,t,x) 
  # approximation of u_t+0.5 u_xx+f(u)=0, u(T,x)=g(x)
  # nonlinearity f and terminal condition g are global
  u = 0
  d = length(x)
  for i = 1 to m^n do
     u = u + g( x+sqrt(T-t)*stdnormal(d) )/m^n
  end for
  for l = 0 to n-1 do
     for i = 1 to m^(n-l) do
        s = uniform(t,T)
        w = x + sqrt(s-t)*stdnormal(d)
        u= u + (T-t)*f(MLP(l,m,s,w))/m^(n-l)
        if l>0 then
           u = u - (T-t)*f(MLP(l-1,m,s,w))/m^(n-l)
        end if
  end for
  return u
end function

# uniform (t,T) returns a sample from the uniform distribution on (t,T)
# stdnormal(d) returns a sample from the d-dimensional standard normal distribution