-\infty
$\begin{array}{r l}{\mathbb{U}[:,,0]}&{{}=~0}\end{array}$
#x-->+\infty
$\#U[:,-1]\ =\ n p.\ e x p\left((1-a)x{-}m a x\ -\ \beta\tau\right)\ <--\$ this would have been the
$\hookrightarrow$ approximated
#from τ we obtain tmarket
tmarket $=\mathrm{~T~}-\mathrm{~(2~~\uptau~/~\upsigma2)}$
$\mathbb{U}\left[:,\mathrm{~\pi{-1}~}\right]=\left(\mathbb{S}{-}\upmu\mathrm{a}\mathbf{x}-\mathrm{~\mathbb{K}~}\mathrm{~\pi{p}~}.\mathrm{exp}\left(-\mathbf{r}\mathrm{~\pi{-}~}\left(\mathbb{T}\mathrm{~-~}\up t{-}\uppi\mathrm{arket}\right)\right)\right)/\mathrm{~\pi{~}\left(\mathbb{K}~}\mathrm{~\pi{-}~}\mathbb{T}\mathrm{~-~}\upmu\left(\mathbb{T}\mathrm{~-~}\mathbf{r}{0}\left(\mathbf{x}-\mathbf{r}{0}\left(\mathbf{x}+\mathbf{r}{-}\right)\right)\right)\right).$
$\hookrightarrow$ $\mathrm{np.exp(\alphax\max\+\\beta\tau)}.$ )

#t $\scriptstyle=O$ U[0,:] $\mathbf{\Sigma}=\mathbf{\Sigma}$ np.maximum(np.exp(((k+1)/2)\x) - np.exp(((k-1)/2)\x), 0)

# 7.a IMPLICIT SCHEME

if chosenmethod $\scriptstyle==$ "implicit":

tridiagonal matrix

delta $\mathbf{\Sigma}=\mathbf{\Sigma}$ dτ / $\mathtt{d x}2$
sub $\mathbf{\Psi}=\mathbf{\Psi}$ -delta $$ np.ones(N - 2)
main $\mathbf{\Sigma}=\mathbf{\Sigma}$ $(1~+~2~$ delta) $$ np.ones(N - 1)
sup $\mathbf{\Sigma}=\mathbf{\Sigma}$ -delta $$ np.ones(N - 2)
#construct the tridiagonal matrix based on the diagonals defined
$\hookrightarrow$ above
A = np.diag(sub, -1) $^+$ np.diag(main) $^+$ np.diag(sup, 1)

#implicit scheme implementation

for j in range(1, $\mathbb{M}{+}1)$ : #rhs i.e.: right hand side rhs $\mathbf{\Sigma}=\mathbf{\Sigma}$ U[j-1, 1:-1].copy() #we take row by row at each iteration #excluding the boundaries (we leave the $\hookrightarrow$ first and #last columns)

#the following two rows correspond to the contributes of the $\hookrightarrow$ boundaries

rhs[0] $+=$ delta $$ U[j,0] #boundary at: x --> -\infty
rhs[-1] $+=$ delta $$ U[j,-1] #bounadry at: and $x\rightarrow+1$ infty
#update the solution of the linear system in the matrix, and
$\hookrightarrow$ we'll use
#it for the next iteration
U[j, 1:-1] $\mathbf{\Sigma}=\mathbf{\Sigma}$ np.linalg.solve(A, rhs) #we solve row by row each
$\hookrightarrow$ unknown variable

# 7.b EXPLICIT SCHEME

if chosenmethod $\scriptstyle==$ "explicit": #definition of the delta delta $\mathbf{\Sigma}=\mathbf{\Sigma}$ dτ / $\mathtt{d x}2$ assert delta $<=\ 0.5$ , f"Stability violated ( $\begin{array}{r}{8=\cdot}\end{array}$ {delta:.3f}>0.5) { $\hookrightarrow$ increase M\ or decrease dx" for i in range(M): #we exclude the first and last column $\mathbb{U}\left[\mathrm{i}{+}1,1{:}-1\right]~=~($ ( delta $$ U[i,0:-2] #take the whole row excluding the last two $\hookrightarrow$ elements $^+$ (1-2\delta) $$ U[i,1:-1] #take the row excluding the first $\hookrightarrow$ and last elements $^+$ delta $$ U[i,2:] #take the whole row starting from the third $\hookrightarrow$ element )

8. TRASFORM BACK TO THE BLACK-SCHOLES MODEL

#trasformation based on the change of variables
$\textsf{V}=$ np.zeroslike(U)
for i in range $(\mathsf{M}{+}1)$ : V[i,:] = K $$ np.exp(α\x + β\τ[i]) $$ U[i,:] #due to how the change of #variables is defined

# 9. PLOTTING SURFACE 3D

Smesh, Tmesh $\mathbf{\Sigma}=\mathbf{\Sigma}$ np.meshgrid(Sgrid, tmarket) #two matrixes composed by
$\hookrightarrow$ the input #arrays repeated
fig $\mathbf{\Sigma}=\mathbf{\Sigma}$ plt.figure() #create the figure object
ax $\mathbf{\Sigma}=\mathbf{\Sigma}$ fig.addsubplot(111, projection $\therefore$ 3d') #adding the 3D subplot to the fig object
ax.plotsurface(Smesh, #x axes Tmesh, #y axes V, #z axes rstride $=5$ , #skip 5 rows for drawing (less heavy to draw) cstride $=5$ , #skip 5 columns for drawing (less heavy to $\hookrightarrow$ draw) cmap $\bullet^{\prime}$ viridis', #heatmap based on the 'viridis scale' alpha $\scriptstyle=0.9$ ) #opacity
ax.setxlabel('S')
ax.setylabel('t')
ax.setzlabel('f(S,t)')