How to implement FEM for a 2D PDE with variable coefficients

Consider the following 2D linear time-dependent PDE with variable coefficients:
begin{equation}label{Bates_equation}
begin{split}
frac{partial u(s,v,tau)}{partialtau} &=frac{1}{2}v s^{2}frac{partial^{2}u(s,v,tau)}{partial s^{2}}+frac{1}{2}sigma^{2}vfrac{partial^{2}u(s,v,tau)}{partial v^{2}}
&quad +rhosigma v sfrac{partial^{2}u(s,v,tau)}{partial spartial v}
&quad +(r-q-lambdaxi)sfrac{partial u(s,v,tau)}{partial s}+kappa(theta-v)frac{partial u(s,v,tau)}{partial v}
&quad-(r+lambda)u(s,v,tau),
end{split}
end{equation}
where $0leq s leq s_{max}=3E$, $0leq vleq 1$, E=100 and some other suitable values.
To solve this PDE using the FEM, I employ the answer given by user21 at FEM differentiation matrices, to construct the FEM differentiation matrices and then obtain the corresponding system of ODEs that must be solved to find the final solution.
I did this as follows:

ClearAll["Global`*"];
Needs["NDSolve`FEM`"]
FiniteElementDerivative[order : {__Integer}, mesh_ElementMesh] /; 
  1 <= Length[order] <= 3 := 
 Block[{dim, nr, vd, sd, mdata, ccoef, pos, dcoef, cdata}, 
  dim = Length[order];
  nr = ToNumericalRegion[mesh];
  vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u}, 
      Table[Unique[X], {dim}]}];
  sd = NDSolve`SolutionData[{"Space"} -> {nr}];
  mdata = InitializePDEMethodData[vd, sd];
  ccoef = ConstantArray[0, dim];
  pos = Flatten[Position[order, 1]];
  ccoef[[pos]] = 1;
  dcoef = ConstantArray[0, dim];
  pos = Flatten[Position[order, 2]];
  dcoef[[pos]] = 1;
  dcoef = DiagonalMatrix[dcoef];
  (*"Pure ConvectionCoefficients" will trigger a warning*)
  Quiet[cdata = InitializePDECoefficients[vd, sd,
     "DiffusionCoefficients" -> {{dcoef}}, 
     "ConvectionCoefficients" -> {{ccoef}}],
   {InitializePDECoefficients::femcscd}];
  DiscretizePDE[cdata, mdata, sd]]

And then write:

m = 32; n = 32; size = m*n;
TT = 1.; r = 0.025;
e = 100.; q = 0; [Sigma] = 0.3; [Kappa] = 1.5; [Theta] = 0.04; 
[Lambda] = 0; [Zeta] = 0; [Rho] = -0.9;
[Rho]1 = Sqrt[1 - [Rho]^2]; [Omega] = r - q - [Lambda]*[Zeta];
xsmin = 0; xsmax = 3 e; ysmin = 0; ysmax = 1.0;

hx = (xsmax - xsmin)/(m - 1); nx = N@Range[xsmin, xsmax, hx]; nx1 = 
 Partition[nx, 1];
hy = (ysmax - ysmin)/(n - 1); ny = N@Range[ysmin, ysmax, hy];
origrid = Flatten[Outer[List, nx, ny], 1];

Idx = SparseArray[{{i_, i_} -> 1.}, {m, m}, 
  0]; xx = (SparseArray@DiagonalMatrix@nx);
Idy = SparseArray[{{i_, i_} -> 1.}, {n, n}, 
  0]; yy = (SparseArray@DiagonalMatrix@ny);
XX = KroneckerProduct[xx, Idy]; Y = KroneckerProduct[Idx, yy];
Id = KroneckerProduct[Idx, Idy];

mesh = ToElementMesh[origrid];
mesh["Wireframe"];

{dXmatFEM, d2XmatFEM} = 
  FiniteElementDerivative[#, mesh]["StiffnessMatrix"] & /@ {{1, 
     0}, {2, 0}};
Map[MatrixPlot, {dXmatFEM, d2XmatFEM}];
{dYmatFEM, d2YmatFEM} = 
  FiniteElementDerivative[#, mesh]["StiffnessMatrix"] & /@ {{0, 
     1}, {0, 2}};
dXYmatFEM = (FiniteElementDerivative[#, mesh][
       "StiffnessMatrix"] & /@ {{1, 1}})[[1]];
Map[MatrixPlot, {dYmatFEM, d2YmatFEM, dXYmatFEM}];

mat = SparseArray@Simplify[(1/2 Y.XX^2).d2XmatFEM
     + (([Rho]*[Sigma])*(XX.Y)).(dXYmatFEM) + (1/
        2 [Sigma]^2 Y).d2YmatFEM
     + ([Omega]*
        XX).(dXmatFEM) + ([Kappa]*([Theta]*Id - 
          Y)).(dYmatFEM) - (r + [Lambda]) Id];
U[t_] = Flatten@Table[Subscript[u, i, j][t], {i, 1, m}, {j, 1, n}];
(*Initial condition*)    
payoff = Flatten@Table[Max[nx[[i]] - e, 0], {i, 1, m}, {j, 1, n}];
initc = Thread[U[0] == payoff];
eqns = Thread[D[U[t], t] == mat.U[t]];

(*Boundaries*)
sf = 0; bc1 = Table[Subscript[u, 1, j][t] == 0., {j, 1, n}];
Table[bc11[l] = Map[D[#, t] + sf # &, bc1[[l]]], {l, 1, Length[bc1]}];
sf = 0; bc2 = 
 Table[Subscript[u, m, j][t] == xsmax*Exp[-q t] - e*Exp[-r t], {j, 1, 
   n}];
Table[bc22[l] = Map[D[#, t] + sf # &, bc2[[l]]], {l, 1, Length[bc2]}];

sf = 0; bc4 = 
 Table[Last[
    NDSolve`FiniteDifferenceDerivative[1, ny, 
     Take[U[t], {(k - 1)*n + 1, (k) n}]],
    "DifferenceOrder" -> 2] == 0, {k, 2, m - 1}];
Chop@Table[
   bc44[l] = Map[D[#, t] + sf # &, bc4[[l]]], {l, 1, Length[bc4]}];

Table[eqns[[i]] = bc11[i], {i, 1, Length[bc1]}];
Table[eqns[[size - i + 1]] = bc22[Length[bc2] - i + 1], {i, 
   Length[bc2], 1, -1}];
Table[eqns[[(i + 1) n]] = bc44[i], {i, 1, Length[bc4]}];

For[k = 2, k <= m - 1, k++,
  bo1 = Drop[Take[U[t], {(k - 1)*n + 1, (k) n}], -1];
  eqns[[(k*n)]] = (Chop@Simplify[
      Last@Table[eqns[[(k*n)]] = eqns[[(k*n)]] /. (D[bo1[[i]], t] ->
            Last@eqns[[((k - 1)*n + 1) + (i - 1)]]), {i, 1, 
         Length@bo1}]
      ]);
  coef1 = Normal@CoefficientArrays[eqns[[k*n]], U[t]];
  coef2 = Coefficient[First[coef1], D[U[t][[k*n]], t]];
  eqns[[k*n]] = 
   D[U[t][[(k*n)]], t] == Chop@(-(1/coef2) Last@coef1).U[t];
  ];

vec0 = SparseArray[{i_} -> 0, size];
mat01 = Table[vec0, {i, 1, n}];
mat02 = Table[-Last@CoefficientArrays[eqns[[i]], U[t]], {i, 
    n + 1, (m - 1) n}];
mat03 = ArrayFlatten[{{mat01}, {mat02}, {mat01}}];
mat03 // MatrixPlot;
vec1 = Last@eqns[[m*n]];
vec2 = SparseArray[{i_} -> 0, ((m - 1)*n)];
vec3 = SparseArray[{i_} -> vec1, n];
vec4[t_?NumericQ] = N@Join[vec2, vec3];
mat03 = SparseArray[Chop@mat03];

(*SOLVING SYSTEM OF ODES*)
Monitor[lines = 
    NDSolve[{D[v[t], t] == mat03.v[t] + vec4[t], 
      v[0] == initc[[All, 2]]}, v[t], {t, 0, TT},
     AccuracyGoal -> 5, PrecisionGoal -> 5,
     Method -> {"FixedStep", "StepSize" -> .0001, 
       Method -> "ExplicitEuler"},
     EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], 
   monitor]; // AbsoluteTiming

s = v[t] /. lines[[1]]; tss = s[[0]][[3]][[1]];
Print["The number of total time steps = ", Length@tss];

s /. {t -> 0};
set0 = Table[Flatten@{origrid[[i]], %[[i]]}, {i, 1, size}];
ListPlot3D[set0, AxesLabel -> {"s", "v", "u"}, PlotRange -> All, 
 ImageSize -> 400]

sol0 = s[[0]][[4]]; sol1 = s[[0]][[4]][[Length@sol0]];
sol2 = sol1[[1]];
list1 = Table[Flatten@{origrid[[i]], sol2[[i]]}, {i, 1, size}];

T12 = Map[Last, list1];
set1 = Table[Flatten@{origrid[[i]], T12[[i]]}, {i, 1, size}];
ListPlot3D[set1, AxesLabel -> {"s", "v", "u"}, PlotRange -> All, 
 ImageSize -> 400]

v0 = 0.04;
f = Interpolation@list1;
Print["The error is = ", ScientificForm@Abs[8.894869 - f[100, v0]]];

Unfortunately, this does not converge to the solution as the last line showing by not approaching zero when I increase the number of nodes, $m$ and $n$.

Can anyone give some hints how I can impose the variable coefficients in the final system matrix? Maybe I am making a mistake at that part! Or the procedure for constructing the matrix in the system of discretized ODEs is different in FEM methodology?

The boundary conditions are:
begin{equation}label{CH3Eq4}
begin{split}
&u(s,v,tau)simeq 0,qquad srightarrow0,
&u(s,v,tau)simeq s_{text{max}}exp{(-qtau)}-Eexp{(-rtau)},qquad srightarrow+infty,
&frac{partial u(s,v,tau)}{partial v}simeq 0,qquad vrightarrow+infty.
end{split}
end{equation}

A sample numerical solution for PDE must look like the following: