Let that we want to numerically solve the following PDE \begin{equation}\label{sde} -r V(S,t)+r S \frac{\partial V(S,t)}{\partial S}+0.5 S^2 \text{sigma}^2 \frac{\partial ^2V(S,t)}{\partial S^2}+\frac{\partial V(S,t)}{\partial t}=0, \end{equation} using the initial and boundary conditions: $V(S,T)=\max (S-K,0),V(1000,t)=1000-K,V(0,t)=0$. So, it would be obvious to write and solve it simply via the following

 ClearAll["Global`*"]
 S0 = 10; T = 1.; sigma = 0.3; r = 0.05; K = 10;
 sol = NDSolve[{D[V[S, t], t] + r*S*D[V[S, t], S]
      + 0.5*sigma^2*S^2*D[V[S, t], {S, 2}] - r*V[S, t] == 0,
    V[S, T] == Max[S - K, 0],
    V[1000, t] == 1000 - K, V[0, t] == 0}, V, {S, 0.1, 1000}, {t, 0, T}]
 Plot3D[V[S, t] /. sol, {S, 0.1, 25}, {t, 0, 1}]

According to "tutorial/NDSolveMethodOfLines", we can use the differentiation scheme by defining the finite difference in matrix notations and solve such PDEs, e.g. as follows:

 S0 = 10; T = 1.; sigma = 0.3; r = 0.05; K = 10;
 a = 0; b = 1.; Subscript[t, 0] = 0; T = 1.; xlen = 1/5; ylen = 1/5;
 xgrid = N[Range[a, b, xlen]]; ygrid = N[Range[Subscript[t, 0], T, ylen]];
 xyGrid =  Table[{xgrid[[i]], ygrid[[j]]}, {i, 1, Length[xgrid]}, {j, 1, Length[ygrid]}];
 dm01 = NDSolve`FiniteDifferenceDerivative[{0, 1}, {xgrid, ygrid}]["DifferentiationMatrix"];
 dm20 = NDSolve`FiniteDifferenceDerivative[{2, 0}, {xgrid, ygrid}]["DifferentiationMatrix"];
 dm10 = NDSolve`FiniteDifferenceDerivative[{1, 0}, {xgrid, ygrid}]["DifferentiationMatrix"];
 wSet = Table[w[j], {j, 1, Dimensions[dm01][[1]]}];

where $wSet$ is the set of unknowns corresponding to the mesh that must be found via solving system of equations. So, if we want to build the numerical procedure, I think we should write as follows:

 Simplify[dm01.wSet + 1/2 sigma^2*S^2*dm20.wSet + r*S*dm10.wSet]

and then solve the resulting system. However, here there are two problems, the first is that how to incorporate $S$ which is related to time in the above discretized systems nd finally how to incorporate the initial and boundary conditions ( $V[S, T] == Max[S - K, 0]$, $V[1000, t] == 1000 - K$, $V[0, t] == 0$)?

I will be thankful if anyone could solve the PDE with the above differentiation scheme or leave me some tips.