Yes. It does the trick. But I want to further achieve the following goals:

1. Show more wave functions, such as 6.
2. Change the displayed style to a more elegant way as shown in the following notebook.

What I mean is, to display these figures in the following way, but at least 6 of them should be displayed:

On the other hand, the notebook represented here has the following problem: It doesn't make sense that the energy eigenvalues, i.e., En[0], En[1], ..., En[4], are artificially set, not calculated from scratch, but I can't think of the code to calculate them automatically.

I tested the time performance and your above code takes over 20 seconds. Is this normal?

In[150]:= startTime = AbsoluteTime[];
Clear[\[HBar], L, \[CapitalGamma], m, \[CapitalEpsilon], x, A, B, G, 
  H, sols];
parameterRules = {\[HBar] -> 1, L -> 1, \[CapitalGamma] -> 100, 
   m -> 1};
\[Alpha] = Sqrt[2 m (\[CapitalGamma] - \[CapitalEpsilon])]/\[HBar];
k = Sqrt[2 m \[CapitalEpsilon]]/\[HBar];
\[Psi]1[x_] = G*Exp[\[Alpha] x];
\[Psi]2[x_] = A*Sin[k x] + B*Cos[k x];
\[Psi]3[x_] = H*Exp[-\[Alpha] x];
l2Norm = 
  Assuming[
   m (\[CapitalGamma] - \[CapitalEpsilon]) > 0 && \[HBar] > 0 && 
    L > 0, Simplify[
    Sqrt[Integrate[\[Psi]1[x]^2, {x, -Infinity, -L/2}] + 
      Integrate[\[Psi]2[x]^2, {x, -L/2, L/2}] + 
      Integrate[\[Psi]3[x]^2, {x, L/2, Infinity}]]]];
regularityConditions = {\[Psi]1[-L/2] == \[Psi]2[-L/2], \[Psi]1'[-L/
      2] == \[Psi]2'[-L/2], \[Psi]2[L/2] == \[Psi]3[L/2], \[Psi]2'[
     L/2] == \[Psi]3'[L/2]};
d = Simplify[
   Det[CoefficientArrays[regularityConditions, {A, B, G, H}][[2]]]];
energyValues[params_] := 
 SolveValues[
  d == 0 && 0 < \[CapitalEpsilon] < \[CapitalGamma] /. 
   params, \[CapitalEpsilon], Reals]
sols[params_] := 
  sols[params] = 
   Table[1/l2Norm*
         Piecewise[{{\[Psi]1[x], 
            x <= -L/2}, {\[Psi]2[x], -L/2 < x <= L/2}}, \[Psi]3[
           x]] /. \[CapitalEpsilon] -> eVal /. params /. 
      Quiet@Solve[
         regularityConditions /. params /. \[CapitalEpsilon] -> 
           eVal][[1]] /. Thread[{A, B, G, H} -> 1], {eVal, 
     energyValues[params]}];
sols[parameterRules] // N // Chop;
Plot[Evaluate[sols[parameterRules]], 
 Evaluate[{x, -L, L} /. 
   parameterRules]]; With[{params = {\[HBar] -> 1, 
    L -> 3, \[CapitalGamma] -> 100, m -> 1}}, 
 Plot[Evaluate[sols[params]], Evaluate[{x, -L, L} /. params]]];
elapsedTime = AbsoluteTime[] - startTime

Out[165]= 22.644842

The way you define ψ[x_] causes an infinite loop of ψ calling itself. Try this:

CostNorm = 1/Sqrt[NIntegrate[Abs[sols]^2, {x, -L, L}]];
normalizedSols = CostNorm*sols;
Plot[normalizedSols, {x, -L, L}, PlotRange -> All]