Remark aside: I often use 'maplemint' and try to understand the answers:

Squaring:=proc(Expr)
...
end proc: 
maplemint(%);

  These names were used as global names, but were not declared:
      _m
  These local variables were never used: 
      E

Strange ... it looks as if the first call defining B is a kind of substitution and the second one evaluates.

So B:=eval(convert(A, list)) gives it, while without eval it does not show it, but knows it. Hm.

Thank you both, I got it now. Neat.

Seems I am too stupid  ... Would you mind to show it again for
Int((sqrt(-x^2+t^2))*t / (exp(t)+1) / Pi^2, t = x .. infinity)
which at x=45 should give ~ 1.1434*10^(-18), now as ODE?

I am afraid that does not work (using t instead of xi to avoid typos):

Int( f(x,t), t=x .. infinity): 'diff(%, x)':
'%' = %;

     /   infinity           \      infinity
     |  /                   |     /
  d  | |                    |    |          d
  -- | |          f(x, t) dt| =  |          -- f(x, t) dt - f(x, x)
  dx | |                    |    |          dx
     |/                     |   /
     \  x                   /     x

So you need that f does not depend on x. 

If you write down your desire and use dsolve you should see it as well.

Note that I corrected my typo, I meant sol1.

I think sol1 can not give the integral (just try without the option 'numeric')
because differentiating Int( f(x,t), t=x .. infinity) w.r.t. x does not give the integrand

Thank you, that what I was looking for.

Concerning the notice: does it mean that in dsolve the integral is used without the intended options?

If I understand you correctly then sol2 sol1 should compute the integral, sol2(45) sol1(45)~ 0.153

But that integral is ~ 1.143 *10^(-18)

rhs(eq2); convert(%, rational);
Int(%, xi=x .. infinity);
eval(%, x=45); evalf(%);

Do you need some complementary value? 0.153 ~ value for x=1

Edited: sol1, not sol2

For values beyond x ~ 36 numerically the denominator is just exp(xi) and
YFD(x) ~ x*(x*BesselK(0,x)+2*BesselK(1,x))/Pi^2 instead of integration:

  Int(sqrt(-x^2+xi^2)*xi/(exp(xi)+1*0), xi = x .. infinity)/Pi^2;
  value(%) assuming 0<x; # gives it

For the first assertion use (the "=" has to be understood 'numerically') 

  exp(xi)+1 = exp(xi); combine(expand(%/exp(xi)));

                           1 + exp(-xi) = 1

But this is almost the definition of DBL_EPSILON

  exp(-xi) <= evalhf(DBL_EPSILON); solve(%, xi);
  
                RealRange(36.0436533891172, infinity)

Yes, MapleCloud works for me - even with my security settings (while the solution pointed to by Markiyan does not [remark aside: they also do use HTML5 canvas data, a technique for fingerprinting a user's browser]). However MapleCloud lives in Amazon's cloud.

@Markiyan Hirnyk 

Due to rounding errors the counter example does not work completely if I
switch to rationals, c <= b+y is not true (Digits=15).

But setting y=Y and using all the other given values one has a polynomial
of degree 2 in Y and solving together with c <= b+Y gives a valid range,
{.234860253652531e-3 < Y, Y < .605803547e-3} in floats for simplicity.

So y = 5/10000 does it.

@Preben Alsholm 

I like it (though knowing almost nothing about integral equations & their numerics), but do not know what happens id switching that to an 'answer'. Minor suggestion: Chebyshev instead of monomials.

@Preben Alsholm 

I am confused now as well, I thought his question is about recursion

So the actual question is to have a numerical approximation for y in 
y(t) = 1-h*Int( T(s) * y(t-s)^4, s = 0 .. t) , 
T = ... having a singularity ~ 1/sqrt(s) in s=0 ?

I was not aware of his earlier post

N = cumulative normal distribution, cN = the complementary of N
If one wants to avoid the usage of imaginary part then one re-writes
using the fact that the erf function and conjugation do commute.

Download MP_invLaplace.mw

One can reduce to phi=1 (see above) and I think a result for 0<x,0<t is

1/4*I*(erf(1/2*(I*2^(1/2)*t+2^(1/2)*t+x)/t^(1/2))*exp(2^(1/2\
)*x+I*t+1/2*I*2^(1/2)*x)+exp(-I*t+2^(1/2)*x-1/2*I*2^(1/2)*x)\
*erf(1/2*(I*2^(1/2)*t-2^(1/2)*t-x)/t^(1/2))-erf(1/2*(I*2^(1/\
2)*t+2^(1/2)*t-x)/t^(1/2))*exp(-1/2*I*(2^(1/2)*x-2*t))-exp(1\
/2*I*(2^(1/2)*x-2*t))*erf(1/2*(I*2^(1/2)*t-2^(1/2)*t+x)/t^(1\
/2))-exp(2^(1/2)*x+I*t+1/2*I*2^(1/2)*x)+exp(-I*t+2^(1/2)*x-1\
/2*I*2^(1/2)*x)-exp(-1/2*I*(2^(1/2)*x-2*t))+exp(1/2*I*(2^(1/\
2)*x-2*t)))*exp(-1/2*2^(1/2)*x);

I do not want to clean that up (dwasons integral or erfi ?) to the Reals.

PS: is that a task from Finance?