# Hybrid ODE simulation with very small parameters

I am trying to solve a Hybrid dynamical system using NDSolve and WhenEvent. I am able to simulate when the parameters are close to 1. However, the practical value of these parameters are close to 0.001. When I put in the actual values, one of two things happens (1) Mathematica crashes, (2) inaccurate simulation.

Parameters:

R = 0.05;G = .04;L = .001;C1 = .001;E1 = 300;Il = 1;us = 0.5;

xd = Inverse[{{R, 1 - us}, {1 - us, -G}}].{E1, 0};
tmax = 2;


Using NDSolve:

    sol = NDSolve[{
x1'[t] == -(R*x1[t] + (1 - u[t])*x2[t] - E1)/L,
x2'[t] == ((1 - u[t])*x1[t] - G*x2[t])/C1,
WhenEvent[x1[t]*xd[[2]] - x2[t]*xd[[1]] > 0.001, u[t] -> 0],
WhenEvent[x1[t]*xd[[2]] - x2[t]*xd[[1]] < -0.001, u[t] -> 1],
x1[0] == 1, x2[0] == 10, u[0] == 0},
{x1, x2, u}, {t, 0, tmax}, DiscreteVariables -> {u}]


Plotting:

Plot[{xd[[2]], x2[t] /. sol, x1[t] /. sol, xd[[2]] - 5,xd[[2]] + 5}, {t, 0, tmax}, PlotPoints -> 10000, PlotRange -> Full]
Plot[u[t] /. sol, {t, 0, tmax}, PlotPoints -> 2000]
Plot[{(x1[t]*xd[[2]] - x2[t]*xd[[1]]) /. sol}, {t, 0, tmax},PlotPoints -> 10000, PlotRange -> Full]


Simulation works well when L=C1=1.

Your code works for me with L = C1 = .001, it's just slow (326 seconds). Looking at the output shows why:
Plot[u[t] /. sol, {t, 0, tmax}, PlotPoints -> 2000]

\$Version