Mechanics for combining likelihood and prior in non-trivial case Bayes rule is simple enough on its face:
$$
pr(B|A) = \frac{pr(A|B)pr(B)}{pr(A)}
$$
If these things are known scalar probabilities, the answer is simple to compute.  
But I'm failing to understand what happens when I'm working with actual probability density functions.  
For example, what if my prior ($pr(B)$) is 
$$
\mathcal{N}\left(\left[\begin{array}{c}{
1\\
1
}\end{array}\right], \left[\begin{array}{cc}
  1 & .5 \\ 
  .5 & 1
 \end{array}\right]
\right)$$
Then I collect some new data that looks like this:
library(MASS)
set.seed(8675309)
N <- 1000

X <- mvrnorm(N, c(0, 0), matrix(c(1, .25, .25, 10), 2))
y <- X %*% c(.5, 2.5) + rnorm(N, sd = 25)

m <- lm(y~X-1)

The likelihood is also normal, with the mean being the coef vector
> coef(m)
       X1        X2 
0.4975087 2.5471237 

and the variance being
> vcov(m) *sigma(m)
           X1         X2
X1 16.4070280 -0.5447927
X2 -0.5447927  1.4990280

How is one to plug this into Bayes rule?  This is what I've been struggling to understand.
Do I simply sample from each of these bivariate normals, multiplying and dividing?  Do I have to worry about the fact that thir elements are correlated?  In other words, am I just doing a hadamard product for multiplying the prior by the likelihood?  And then the same thing in inverse to scale?
A simple simulation makes it seem like that's not the case; these results are weird:
AB <- mvrnorm(N, coef(m), vcov(m) *sigma(m))
B <- mvrnorm(N, c(1,1), matrix(c(1,.5,.5,1),2))
A <- mvrnorm(N, c(0, 0), matrix(c(1, .25, .25, 10), 2))

bhat <- AB*B/A

hist(bhat[,1])
hist(bhat[,2])


And where does MCMC come in?  Is that just a way to figure out where in the distributions I should be sampling from?
Would appreciate an answer that could help me fill in the gaps!
EDIT
I think that what I needed was for someone to tell me "do it over a grid" or "do it over the domain of your data."  And "multiply the densities, point by point, not the values"
Here's a univariate example:
G <- seq(-5, 5, .01)
pB <- dnorm(x = G, mean = -3, sd = .5)
pB <- pB/sum(pB)
plot(G, pB, type = "l")
pAB <- dnorm(x = G, mean = 3, sd = .5)
pAB <- pAB/sum(pAB)
lines(G, pAB, type = "l", col = "red")
pBA <- pB*pAB/(sum(pB*pAB))
lines(G, pBA, type = "l", col = "blue")

But I must be doing something wrong here:  why is the posterior pointier than the prior and the likelihood in this case?  Or is it that the normalization something more than a factor to scale the PDF to sum to 1?

But I must be doing something wrong here:  why is the posterior pointier than the prior and the likelihood in this case?  Or is it that the normalization something more than a factor to scale the PDF to sum to 1?
 A: I think your problem setup is a little troublesome. Let's take a step back and write down a model.
Bayes' Rule says
$$ p(\theta \vert y) \propto p(y \vert \theta) p(\theta)$$
Equivalently
$$\log(p(\theta \vert y)) \propto \log(p(y \vert \theta)) + \log( p(\theta))$$
I'm going to compute things on the log scale because it's easier.
Let's say I have some data that I think comes from a bivariate normal with known covariance matrix (assume it is the identity for simplicity). Our job is to estimate the mean of the bivariate normal from data. Our likelihood is then
$$ y \vert \mu \sim \mathcal{N}(\mu, \Sigma) $$
Now we need a prior for mu.  Let's just use the one you gave
$$ \mu \sim \mathcal{N}\left(\left[\begin{array}{l}{1} \\ {1}\end{array}\right],\left[\begin{array}{cc}{1} & {.5} \\ {.5} & {1}\end{array}\right]\right) $$
Now, we need to evaluate the posterior along a 2d grid of possible values of $\mu$.  To compute the posterior probability (well, not really, just something that is proportional to the probability) we compute
$$\log(p(y \vert \mu)) + \log( p(\mu))$$
But the important part is that we need to compute this over a grid.  That grid will be the elements of $\mu$.  Let's turn to python to do this
import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt
#Data Generating Parameters
true_mean = [-1.5, -1.75]
true_Sigma = np.eye(2)

#Generate data from our data generating process
data = multivariate_normal.rvs(mean = true_mean, cov = true_Sigma, size = 100)

#Our prior
prior = multivariate_normal(mean = np.ones(2), cov = np.array([[1,0.5], [0.5,1]]))

#Generate a grid on which to evaluate the posterior
grid = np.linspace(-2,2,26)
mu_1, mu_2 = np.meshgrid(grid,grid)
#Reshape the grid for ease of computation
MU = np.stack((mu_1, mu_2), axis = -1).reshape(-1,1,2)

#Instantiate the loglikelihood
loglik = prior.logpdf(MU)

fig,ax = plt.subplots(dpi = 120)
ax.axes.set_aspect('equal')
plt.contourf(mu_1, mu_2,loglik.reshape(mu_1.shape))
plt.title('Prior')
plt.show()

LL = []
for m in MU:
    #Evualuate the likelihood of the data over the grid
    LL.append(multivariate_normal.logpdf(data, mean = m[0], cov = np.eye(2)))

LL = np.array(LL)    
#Add the log likelihood to the log prior
loglik+=LL.sum(axis = 1)



fig,ax = plt.subplots(dpi = 120)
ax.axes.set_aspect('equal')
plt.contourf(mu_1, mu_2,loglik.reshape(mu_1.shape))
plt.title('Posterior')
plt.show()

This code produces the following


As you can  see, the prior looks like your prior distribution (as it should).  After seeing data, the posterior has shifted the density towards the true mean of the data generating process.  So it isn't enough to just evaluate bayes rule as you do in your samples.  You have to work with densities, which can get tricky.

And where does MCMC come in? Is that just a way to figure out where in the distributions I should be sampling from?

As you can imagine, the densities become very difficult to deal with with a modest amount of parameters. MCMC and HMC (Hamiltonian Monte Carlo) are ways to generating markov chains with limiting distributions equal to our posterior. There is a lot to say about MCMC and HMC, so I suggest you create another question if you are really interested.
