There is a pretty short proof (usually called the martingale proof), once you established some major theorems. In particular, we assume that we know the Fundamental theorem of asset pricing and some properties of brownian motions.
So, BS start assuming that the stock dynamic is
$$\frac{dS_t}{S_t} = \mu dt + \sigma dB_t$$
Using Ito's lemma, you can show that the solution to this SDE is
$$S_t = S_0e^{(\mu - \sigma^2/2)t + \sigma B_t}$$
i.e. $S_t$ follows the so called geometric brownian motion.
Now, we show that the market defined by $S_t$ and the risk free rate is complete (intuitively, there is only one source of randomness, in this case the brownian motion $B_t$). If the market is complete, it means that everything can be prices by discounting the expectation of the payoff under the risk neutral measure; for the case of a call option, this is
$$V(t, T) = e^{-r(T-t)} E_{Q}[(S_T - K)^+ \mid \mathcal F_t]$$
Note that we are using the risk netrual measure $Q$ instead of $P$, the "real world" probability measure. This means that we don't know (yet) the dynamic of $S_t$ under the $Q$ measure.
How to find the dynamics under $Q$? Well, we start from the observation that under $Q$, the drift of the stock price $S_t$ must be equal to the risk free rate $r$. We will also use Girsanov's theorem, that says that if
$$\frac{dQ}{dP}\mid _{\mathcal F_t} = \mathcal E(B_t) = \exp\left(\alpha B_t - \frac 12 \alpha^2t\right)$$
then $dW_t = dB_t - \alpha dt$ is also a brownian motion. Substututing in the initial SDE and imposing that the drift is $r$, we get $\alpha = \frac{r - \mu}{\sigma}$ and the dynamics of $S_t$ under $Q$ are
$$S_t = S_0 e^{(r - \sigma^2/2)t + \sigma W_t}$$
Notice how $\mu$ is not a part of the equation anymore! And in fact the price of an option is independent on the drift of the stock under the real-world measure $P$, pretty counter-intuitive at first.
Now, it's almost done! We just need to compute an expectation ;)
$$V(t,T) = e^{-r(T-t)} E_{Q}[(S_T - K)^+ \mid \mathcal F_t]$$
and we know the dynamics of $S_T$ under $Q$. If for simplicity we take $t = 0$ (otherwise we need to condition use the fact that $S_t$ is known when we condition wrt to $\mathcal F_t$.. nothing really changes). We just need to compute
$$\begin{align} V(0,T) = e^{-rT} E_{Q}[(S_T - K)1_{S_T > K}] &= \\ e^{-rT} E_{Q}[S_T1_{S_T > K} - K1_{S_T > K}] &= \\ e^{-rT}E_{Q}[S_T1_{S_T > K}] - e^{-rT}KQ(S_T > K) \end{align} $$
Now $Q(S_T > K)$ is trivial to compute;
$$Q(W_T > \log \frac K{S_0} - \frac{r - \sigma^2}{2} T) = N(d_2)$$
because $W_T$ (Under Q) is normally distributed with variance $T$.
For the first term there is a similar trick but I don't remember at the moment.. But you get the idea, now you just have a function of a normal random variable (i.e. $S_T$) and finding the expectation should be easy
Best Answer
You shouldn't start with a PDE. The Black-Scholes-Merton equation is a PDE that lists the relationship between the greeks. You will need to differentiate with the solution of the PDE . You have listed the solution to the PDE for a call option in your first part of the question. You will need to do partial differentiation for the solution. $\frac{\partial^2 C}{\partial S^2}$ is the gamma of the option.
I could derive those for you but there are tons of material online that you can follow. Just remember, the PDE doesn't tell you anything without a boundary condition. The fact that you have listed the PDE but without a boundary condition for the call option isn't technically correct. Once you have a boundary condition, you'll be able to derive a solution. This solution is your first part of the question.