How to improve the solution of the following differential equation?

Consider a PDE
$$
tag 1 frac{partial f_{n}}{partial t}- frac{1}{2t}Efrac{partial f_{n}}{partial E} = 0
$$
I choose the initial condition $f_{n}(E, t_{0}) = e^{-E/T_{0}}$, where $T_{0}$ is some constant (plus the boundary condition $f_{n}(E_{max},t) = 0$ for some large $E_{max}$).

Eq. (1) may be transformed into a trivial ODE by performing a change of variables. For the given initial condition, the $f_{n}(E,t)$ is
$$
tag 2 f_{n}(E,t) = e^{-E/(T_{0}sqrt{t0/t})}
$$
I want to obtain this result by solving the PDE. This is my code: first, I discretize the E domain, then turn the PDE into ODE by discretizing the derivative, and then solve the set of ODEs:

imax = 500;
t0 = 0.1;
(*E grid*)
DEvalue = 0.05;
Emin = 0.01;
tmax = 4;
EStep[k_, DE_] = Emin + DE*k;
(*Discretization of the derivative*)
fnEnDerivative[i_] = 
  If[i != imax, (fn[i + 1][t] - fn[i][t])/DEvalue, -fn[imax][t]/
    DEvalue];
Tt[t_]=0.84/Sqrt[t];
(*Table with equations, initial conditions and functions*)
InitialConditionTable = 
  Join[{Tg[t0] == Tt[t0]}, 
   Table[fn[i][t0] == Exp[-(EStep[i, DEvalue]/Tt[t0])], {i, 0, imax, 
     1}]];
functionsTable = Table[fn[i], {i, 0, imax, 1}];
EquationsTable = 
  Table[fn[i]'[t] -
     1/(2*t)*EStep[i, DEvalue]*
      fnEnDerivative[i]== 0, {i, 0, 
    imax, 1}];
sol = NDSolve[{EquationsTable, InitialConditionTable}, 
    functionsTable, {t, t0, tmax}, 
    Method -> {"EquationSimplification" -> "Solve"}][[1]];

After obtaining the solution, I compare the behavior of the exact solution $(2)$ and the numerical solution. I found a discrepancy due to low precision:

fnv[t_] := 
     Table[{EStep[i, DEvalue], (fn[i] /. sol)[t]}, {i, 1, imax, 1}]
    BoltzmannDistrDiscrete[t_] := 
     Table[{EStep[i, DEvalue], Exp[-EStep[i, DEvalue]/Tt[t]]}, {i, 1, 
       imax, 1}]
    ListLogPlot[{fnv[0.3], BoltzmannDistrDiscrete[0.3]}, 
     PlotRange -> {{0.003, 25}, All}, Joined -> True, Frame -> True, 
     ImageSize -> Large]

Improving Accuracy/Precision Goals does not help. For example, with PrecisionGoal ->30 I get this:

Could you please help me with how to resolve this problem?

imax = 500;
t0 = 10/11;
(*E grid*)
DE = 1/20;
Emin = 1/100;
grid = Table[Emin + DE*k, {k, 0, imax}];

fd = NDSolve`FiniteDifferenceDerivative[Derivative[1], grid, 
  "DifferenceOrder" -> 8]; m = fd["DifferentiationMatrix"];
var = Array[f, {Length[grid]}];
vart = Table[f[i][t], {i, Length[grid]}]; fx = m.vart;

eqs = Table[
   f[i]'[t] - 1/(2 t) grid[[i]] fx[[i]] == 0, {i, Length[grid]}];

Tt[t_] = 84/100/Sqrt[t];
(*Table with equations,initial conditions and functions*)
InitialConditionTable = 
  Table[f[i][t0] == Exp[-(grid[[i]]/Tt[t0])], {i, Length[grid]}];

tmax = 4/10; sol = 
 NDSolve[{eqs, InitialConditionTable}, var, {t, t0, tmax}, 
   Method -> {"EquationSimplification" -> "Solve"}][[1]];
fnv[t_] := Table[{grid[[i]], (f[i] /. sol)[t]}, {i, Length[grid]}]
BoltzmannDistrDiscrete[t_] := 
 Table[{grid[[i]], Exp[-grid[[i]]/Tt[t]]}, {i, Length[grid]}]
ListLogPlot[{fnv[tmax], BoltzmannDistrDiscrete[tmax]}, 
 PlotRange -> {{.003, 25}, All}, Joined -> True, Frame -> True, 
 ImageSize -> Large]

tmax = 4/10; sol = 
 NDSolve[{eqs, InitialConditionTable}, var, {t, t0, tmax}, 
   Method -> {"EquationSimplification" -> "Solve"}, 
   WorkingPrecision -> 30][[1]];