Well, the problem is in $P$ after all.
I'll keep the previous answer as it also works for the directed case (which is NPC, as answered on the other question), and shows it is $FPT$ with respect to $l$.
In the undirected case, it is solvable, deterministically via minimum cost flow (this might not work on the scales you are referring to in the question, but its better than exponential algorithm.
The following procedure will decide whether some edge $e=(u,v)\in E$ should be a part of the output graph. In order to answer the original problem just loop over all edges.
In order to create the flow network, do as follows:
Step 1:
Expand $e$ to have a vertex $x_e$ and replace $e$ with the edges $(u,x_e),$$(x_e,u),(v,x_e),(x_e,v)$ (they are directed as a part of the flow network), set their cost to 0.
Step 2: replace every vertex $t$, except for $x_e$ by two vertices $t^-$ and $t^+$, and add an edge $(t^-,t^+)$. Set the cost of these edges to 1.
Step 3: Replace of every edge $\{a,b\}\in E$ with the edges $(a^+,b^-),(b^+,a^-)$. Set the cost of these edges to 0.
Step 4: Add a new vertex $y_e$ and add the edges $(s,y_e),(t,y_e)$ with cost 0.
Step 5: set all capacities to 1.
Now run the min cost flow algorithm, searching for a flow of value 2 from $x_e$ to $y_e$.
Analysis:
- Every 2-valued flow from $x_e$ to $y_e$ is a union of a path $x_e\leadsto s\to y_e$ and a path $x_e\leadsto t\to y_e$.
- The paths are disjoint, since for every vertex $t$ there's only 1 capacity in the $(t^-,t^+)$ arc.
- The returned paths are the two paths whose sum of distances is minimal, and that's also the cost of the found flow. This allows us to add $e$ to the output graph or delete otherwise.