Joint expectations in Python or R $X$ and $Y$ are Bivariate Normally distributed as $N(0,1)$ with a correlation coefficient of $p$.
Is there any way to find the expectation of $f(x)*g(y)$ - Let $f(x)$ this be a function of $x$ which is Integrable as well as Differentiable
$E[f(x)*g(y)]$? Where $g(y) = 1$, if $y>c$. $g(y) = 0$, otherwise
Is there any way to do this in Python? Or in some other Language like R?
I am aware of the multivariate_normal() function under scipy.stats. 
import scipy.stats as st

a = st.multivariate_normal()

Is there any way we can use this combined with some other method?
It will be great if someone can help me with an analytical solution also. For example we can take $$ f(x) = exp(x)$$
 A: You can rather easily approximate the expected value of $f(X)g(Y)$ using simulation. For instance, here is the R code to simulate the expected value of $X^2e^Y$ with $\rho = 0.2$ using 1 million samples:
# Parameters for calculation
f <- function(x) x^2
g <- function(y) exp(y)
rho <- 0.2

# Approximate expectation
library(mvtnorm)
set.seed(144)
simulated <- rmvnorm(1e6, c(0, 0), rbind(c(1, rho), c(rho, 1)))
mean(f(simulated[,1]) * g(simulated[,2]))
# [1] 1.713411

Similar code can approximate the expectation in python:
import numpy as np
np.random.seed(144)
simulated = np.random.multivariate_normal([0, 0], [[1, 0.2], [0.2, 1]], 1000000)
print(np.average(simulated[:,0] ** 2 * np.exp(simulated[:,1])))
# 1.70735583203

Analytically, the expectation is given by the following expression, based on the bivariate normal pdf:
$$
\int_{-\infty}^\infty \int_{-\infty}^\infty f(x)g(y)\frac{1}{2\pi\sqrt{1-\rho^2}}\exp\bigg[-\frac{x^2-2\rho xy + y^2}{2(1-\rho^2)}\bigg] dxdy
$$
Whether this has an analytical solution depends on the functions $f$ and $g$.
