One can simplify the integrals over the complex error function into just one, which speeds up.

I did not use acer's error handling / code (I do not quire understand it), but a cut off.

A proper way for the nasty error  function would be through https://en.wikipedia.org/wiki/Faddeeva_function because of the damping by exp(-kappa^2). There is code for it (through the literature on the German Wikipedia site), but not in Maple.

Edit 13.09.22: the cut off = 6 results from roughly estimating an (absolute) integration error to be less than 1e-15.

Download MP_234736.mw

For example:

Download MP_234721.mw

Essentially you need Int(x^(1/10)*exp(x),x=0..2) = 6.47950080882301

Int(x^(1/10)*exp(x),x);
IntegrationTools:-Parts(%, exp(x));
value(%):
V:=simplify(%) assuming 0<x;

Plotting "shows" it is continous [Re(V), Im(V)]: plot(%, x=-0.1..2.1);

limit(V, x=2, left)-limit(V, x=0, right);  
expand(%);
evalf(%);

"proves" the value
-1/10*(-1)^(9/10)*GAMMA(1/10,-2)+2^(1/10)*exp(2)+1/10*Pi*(-1)^(9/10)/sin(1/10*Pi)/GAMMA(9/10)
is correct

Have not tried for your final integral

You want to write proc(x,y,t), else y is local (as you did) and not known outside

Cant you feed the solution into the ode and then using MultiSeries:-series ?

Download MP_234443_odetest_series.mw

You have many terms, 196*109. Here is a suggestion: write the inner Sum as function, evaluate the outer additions and then feed the function with its numerical values.

I changed your variable l, because it is bad to read it and reduced addition to 20*10 (well, I lost patients). Also note that you may have to do some more work if you want to care for rounding errors of 20000 terms of quite different magnitude and signs.

Hope it helps.

MP_234318.mw

Edit: I need roughly 200 sec to get 5.76900063193468*10^29 with an old Win PC, summing up over your original range.

Likewise you can look into the Digital Library of Mathematical Functions https://dlmf.nist.gov/,

https://dlmf.nist.gov/19.6

I would do it in a paper & pencil style (it is a mess how this board displays it)

Download MP_234277_CubicSolve.mws

MP_234277_CubicSolve.pdf

As far as I see you have a linear system 4 x 4. The according matrix has rank <= 3. Hence the solution is not a point

MP234160_QuestionNo1.mw

For your example the following would do:

sin(x)/x; convert(%, hypergeom, include=sin); convert(%, BesselJ);

or even shorter

sin(x)/x; convert(%, BesselJ, include=sin);

@vv 

Sum(Beta(k,1/2)/(2*k+1)^2,k=1..infinity);
convert(%, MeijerG);
#convert(%, hypergeom); # needed if evalf is used
convert(%, Int): value(%);
                        infinity
                         -----
                          \      Beta(1/2, k)
                           )     ------------
                          /                2
                         -----    (2 k + 1)
                         k = 1

          1/2
    1/4 Pi    MeijerG([[-1/2, 0, 0], []], [[0], [-3/2, -3/2]], -1)

                            4 - 4 Catalan

You may try assuming, not assume. Or use "simplify(%) assuming ..." for the result.

BTW: Maple 2021 returns /(A^2-B^2)^(1/2)*EllipticK(B*2^(1/2)*(-1/(A^2-B^2))^(1/2))

I posted an answer at https://www.mapleprimes.com/questions/233780-Plot-Absolute-Advantage-Of-Option-A   but that thread disappeared (was it deleted?)

VT:=V_TEV-V_P:
indets(%, symbol);
                                {E, T}

So your VT depends on E = income and T = time and the following shows the values:

plot3d(VT, E=10^5 ... 10^6, T = 1 ..5, axes=boxed);

plot3d(VT/E, E=10^5 ... 10^6, T = 1 ..5, axes=boxed); shows you the values relative to income

Remark: your r seems to stand for rates, 0.1 = 10% is somewhat high, perhaps you want to see it for 1% or 3%

You may consider to use Excel and - if really needed - connect Excel and Maple