Natural boundary conditions vanish in your variational formulation. The following works as expected.
problem Laplace(u, v, solver=CG) = int3d(Th)(Grad(u)' * Grad(v)) + on(1, u=20.0) + on(5, u=-13.1);