From Mathematica symbolic PDE document.
Solve for \(V(S,t)\) where \(V\) is the price of the option as a function of stock price \(S\) and time \(t\). \(r\) is the risk-free interest rate, and \(\sigma \) is the volatility of the stock. \[ \frac {\partial V}{\partial t} + \frac {1}{2} \sigma ^2 S^2 \frac {\partial ^2 V}{\partial S^2} = r V - r S \frac {\partial V}{\partial S} \] With boundary condition \( V(S,T) = \max \{ S-k,0 \}\)
Reference https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_equation See the European call version at bottom of the page.
Mathematica ✓
ClearAll["Global`*"]; pde = D[u[x, t], t] == (1*sigma^2*D[u[x, t], {x, 2}])/2; ic = u[x, 0] == k*Exp[x - 1]*HeavisideTheta[x]; sol = AbsoluteTiming[TimeConstrained[DSolve[{pde, ic}, u[x, t], {x, t}, Assumptions -> k > 0], 60*10]];
\[\left \{\left \{u(x,t)\to \frac {1}{2} k e^{\frac {\sigma ^2 t}{2}+x-1} \left (\text {Erf}\left (\frac {\sigma ^2 t+x}{\sqrt {2} \sqrt {t} | \sigma | }\right )+1\right )\right \}\right \}\]
Maple ✓
restart; pde := diff(u(x,t),t) = 1/2*sigma^2*diff(u(x,t),x$2); ic := u(x, 0) = k*exp(x - 1)*Heaviside(x); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([pde,ic],u(x,t)) assuming k>0),output='realtime'));
\[u \left ( x,t \right ) =k{{\rm e}^{-1}} \left ( -i{\it invfourier} \left ( {\frac {1}{s+i}{{\rm e}^{-{\frac {{s}^{2}{\sigma }^{2}t}{2}}}}},s,x \right ) +{\it invfourier} \left ( {{\rm e}^{-{\frac {{s}^{2}{\sigma }^{2}t}{2}}}}{\it fourier} \left ( {{\rm e}^{x}},x,s \right ) ,s,x \right ) \right ) \]
____________________________________________________________________________________