@dharr 

I have implemented the solution in Mathematica using a modular approach with NumericQ guards to ensure robust numerical evaluation. The integration is split into the diagonal 'wings' ($\lambda_1, \lambda_2$) and the horizontal segment ($\lambda_3$)
approxW[x_?NumericQ, t_?NumericQ] :=
Module[{temp1, temp2},
  temp1 = NIntegrate[
    Re[v[λ1[r], x, t]*dλ1[r] - v[λ2[r], x, t]*dλ2[r]], {r, 0, Infinity},
    AccuracyGoal -> 8, PrecisionGoal -> 8, MaxRecursion -> 25];
  temp2 = NIntegrate[
    Re[v[λ3[s], x, t]*dλ3[s]], {s, -l, l},
    AccuracyGoal -> 8, PrecisionGoal -> 8, MaxRecursion -> 25];
  temp1 + temp2]

approxW[0.01, 0.1] = 0.375077, which is different from my solution

@dharr 

The Fokas Method represents the solution as an integral in the complex $\lambda$-plane. To ensure the integral converges and is well-defined, we must stay within the $D^+$ domain, where the exponential term $e^{-d\lambda^2 t}$ decays rapidly.

Avoiding Poles: The boundary conditions (Robin type) introduce "poles" (points where the denominator becomes zero) along the real axis. By lifting the horizontal contour $\lambda_3$ by $+i$ (the $s+i$ shift you used), you effectively "fly over" these singularities, making the numerical integration smooth and reliable.

Analyticity: Because the function is analytic within the region enclosed by your chosen contours ($\lambda_1, \lambda_2, \lambda_3$), the Cauchy Integral Theorem ensures that this path provides the exact solution for the heat equation part of the Hopf-Cole transform.

I have solved this problem with the Fokas method and with the Fourier Series expansions, and as you can observe, the results are similar but not the same. For this reason I asked how to create  more accurate plots

Download ex2.mw

Download Fourier_ex.2_.mw

@dharr @acer Thanks for your assistance. Should I increase the tolerance for more accurate plots?

@dharr First of all, thank you very much for your help. Numerically, I validated my research, although when I tried to create the plots, Maple ran for too much time. Is there a trick to creating the plots faster? I have highlighted the code blocks in red.

Download third_try.mw

@dharr I have tried what you suggested, but as you can see, approx_w is not numeric. It should be 0.375077

Download second_try.mw

@acer 

Yes exactly! I tried running it again after fixing k*(L-x), but it still runs indefinitely.

Download resources.mw

@sand15 You are right the integrand must be integrand1 := Re(eval(V, k = k1(r))*dk1(r) - eval(V, k = k2(r))*dk2(r)):

@dharr Thank you!! It was difficult to find what to do hear. There were no other references for Fokas Method on Maple!

@Rouben Rostamian, I created a function, so the solution is given faster now, but when I try to replicate the two panels from the textbook, as you can see they are a bit different

Download fokas_method_new.mw

@dharr Thank you for your answer and the time you spent looking at my question. I will try to find a faster way. 

@dharr oh!! Thank you this info is very usefyl.

Download Plotting_contour_integrals.mw

@mmcdara  I played a bit with the code you gave me, and now I have the desired results. Thanks

@Rouben Rostamian  I thought about that, and sure, it's handy. However, I would like to know if I can do the same with Maple. I searched on Google for a tutorial, but I did not find anything.

@mmcdara, sorry for my English!! This is exactly what I was trying to do