I follow a recipe given by Preben Alsholm, see the fork given at
http://www.mapleprimes.com/questions/204941-Dsolvenumeric-And-Numerical-Functions

The following is *very* slow, but it works. So for plotting one may wish to
consider the option 'numpoints=5' or similar.

f:= x -> evalf( Int((sqrt(-x^2+xi^2))*xi / (exp(xi)+1),xi = x .. infinity, method = _d01amc)/Pi^2 );
F:= proc(x) if not type(x,numeric) then 'procname(_passed)' else f(x) end if; end proc;

ODE:= diff(Y(x), x)+10^5*(Y(x)^2-F(x)^2)/x^2;
IC := Y(1) = F(1);

sol:=dsolve({ODE, IC}, numeric, known=F, method='gear');
S:= x -> rhs(sol(x)[2]);

Using Digits = 15 I get

S(1.1) = 0.148400199613981
S(15)  = 0.000191332950036288

And made the experience, that one should not interrupt calculations, since
that may give false values in calling fct values 'near by'

Asymptotically for x beyond ~ 36 one may speed it up. I guess.

I think your recursion writes as

y(k+1,t) = 1-h*Int(JacobiTheta3(0,exp(-1/s))/(Pi*s)^(1/2)*y(k,t-s)^4,s = 0 .. t)

I also doubt you can find the (numerical) limit w.r.t. k (by the way: for which t?)

Pi, not pi is what you should use

Option Explicit

Function W(ByVal z As Double) As Double
Dim x As Double, x_new As Double, NewtonStep As Double
Dim i As Integer

If (z < -1# / Exp(1#)) Then
  W = 1# / Sin(0#)
  Exit Function
End If

i = 0
If (z <= 2#) Then
  x = 0.5 * z
Else
  x = Log(z)
End If
For i = 0 To 8
  If (z <= 2#) Then
    NewtonStep = x - (x * Exp(x) - z) / (1# + x) / Exp(x)
  Else
    NewtonStep = x * (1# + Log(z / x)) / (1# + x)
  End If
  x_new = NewtonStep
  If (CDbl(x) = CDbl(x_new)) Then
    Exit For
  End If
  x = x_new
  i = i + 1
Next
W = x
End Function


Sub tst_W()
Dim x, y
x = 2
Debug.Print W(x)

y = 12.345
x = y * Exp(y)
Debug.Print x, W(x), y

y = -0.9999995
x = y * Exp(y)
Debug.Print x, W(x), y
End Sub

To finish it symbolically I consider acer's EEb as a Laplace transform:

'inttrans[laplace](-1/5*Sin(w)^3/w^(3/5),w,t)': 
'%'=combine(convert(%, Int));
                               3
                         Sin(w)
  inttrans[laplace](-1/5 -------, w, t) =
                           3/5
                          w

           infinity
          /                    3
         |               Sin(w)  exp(-w t)
         |          -1/5 ----------------- dw
         |                      3/5
        /                      w
          0

-1/5*sin(w)^3/w^(3/5):
inttrans[laplace](%,w,t): 
limit(%, t=0, right): simplify(%); evalf[17](%);

                          1/2     3 Pi    (3/5)
                      Pi 2    cos(----) (3      - 9)
                                   10
                 1/30 ------------------------------
                                          1/2 1/2
                         GAMMA(3/5) (5 + 5   )


                         -0.15356212695352348

I rarely use that feature in Maple. Just staring at your picture and its symmetries I perhaps would generate the Gaussian - and 2 indicator functions on the intervals (being 0 or 1), one for the tails and another around zero.

Multiplication should give what you want (but which extra conditions?)

Indicators are easy: generate 2 random floats and use 'piecewise'

Does that help?

seqLimit:=proc(expr,n,ord:=2)
  local x;
  subs(n=floor(x), expr);
  MultiSeries:-asympt(%, x, ord);
  convert(%, polynom);
  simplify(%) assuming x::posint; # lprint(%);
  limit(%, x=infinity);           # just in case
end proc;

  sin((n^2+1)/n*Pi);
  seqLimit(%, n);
 
                                  0

  x+sin((n^2+1)/n*Pi);
  seqLimit(%, n, 3);

                                  x

  sin(sqrt(n^2+n)*Pi)^2;
  seqLimit(%, n);

                                  1

So you ask to solve

[2*f1*r^4, 2*f2*r^2]; L:=simplify(%, size);

  [(-a^2-3*a)*r^(a+2)+(-2*n-2)*r^(n+2)-4*q*r*phi[0]+2*r^(c+2)+1 = 0,   (2*a^2+a)*r^(a+2)+(n^2+n)*r^(n+2)-2*p*r*phi[0]-r^(c+2)*c+1 = 0]

for q, r, given a, n, c, phi[0], p.

subs(n=n-2, a=a-2, c=c -2, L); eliminate(%, q); simplify(%, size);

This gives: q ist is linear in the others and you have to solve

(2*a^2-7*a+6)*r^a+(n^2-3*n+2)*r^n+(-c+2)*r^c-2*p*r*phi0+1

for r. I doubt there is a symbolic solution, even for Reals, thus you want to use numerical solutions (if they exist).

That says that this construct is not supported. Perhaps I would try 'arrays'.

Note the Upper Case for Sum (and you do not use a range?)

2*(Sum(A*cos(omega*t)^2+cos(omega*t)*B*sin(omega*t)-cos(omega*t)*ymeas(t), t));
combine(%);
op(1, %);
[op(%)];
map('q -> Sum(q,t)', %);
convert(%, `+`);

After Carl's correction: the inner integral PT1-PT2 has a value, but not each of the summands if they would stand alone and I think the value is Inner := Pi*p/(b^2+2*abs(b*m)+m^2+p^2)/abs(b).

Then the final integral is rapidly found by Maple, now with dependencies on signs for m and b which makes it look a bit complicated.

For the case 0 < m, 0 < b it is quite simple, -1/4*c*Pi^3/b^2/(2*b+m)^2.

You may wish to verify numerically by choosing values for the parameters, "Inner" is enough at first step.

Edited: have not looked deeper into b=0.

Using Carl's "non-sloppy form" one sees that the task is linear in b,

H=Int((`@@`(D,2)(f)(Y)*Int((-X*b+1)*Int(D(f)(x)^2,x = 0 .. X),X = 1 .. Y)
  +D(f)(Y)*(-Y*b+1)*Int(D(f)(x)^2,x = 0 .. Y))*f(Y),Y = 0 .. 1), f = ms

So 2 results for 2 values of b should do, using his suggestion.


My route was different: I transformed the task to the cube [0 .. 1]^3.

It turns out that f, f', f'' have inputs in [0 .. 1] after that (on the now
combined variables), plotting gives nice shapes in that range.

I used chebpade to approximate f by a polynomial (real + imag*I). For which
one can expect that H has a solution (either wait or help Maple for inner
integrals, degrees may be high and coefficients=rational become large).

Finally I get 

Sol := b -> 26.0856584458258 + .616281090328306e-7*I 
  -(21.9402269803097 + .514632049063567e-7*I)*b

Which re-decovers linearity in b.

Sol(7/10): evalf(%);
                                                    -7
             10.7274995596090 + 0.256038655983809 10   I


PS: it would be better to approximate f' and f'' separately.

By "stomach" says that is something like a bi- or tri-variate cumulative normal distribution. And there is no solution in terms of 'other' functions: Integral(exp * erf) ~ bivariate and you have some limiting case, no?

for that very example it is easy to see that there is no real solution:

Your first equation in your system S is a polynomial in x3 only. Now look up 'sturm' in the help to see that the following answer confirms what I have asserted:

sturm(S[1],x3,-infinity,infinity);

                                  0

Of course not as elegant as Rouben Rostamian's neat solution

As it was mentioned consider the general case, instead only k=5. Again split in the zero of the integrand, call it b. Then b^k = 1/2 and for the integral found by Maple in terms hypergeometric 2F1 one has to prove

hypergeom([1/k, -1/k],[(k+1)/k],1/2) =   2^(-(k+1)/k)+2^(-(k-1)/k)*GAMMA((k+1)/k)^2/GAMMA((2+k)/k)

Now rewrite that 2F1 as integral and then recognize that as an incomplete Beta function value.

Using http://dlmf.nist.gov/8.17, identity 8.17.21 then "proves" it.

Edited: please find a worksheet appened

MP_integral_with_absolute_and_hypergeom_beta_(2).mws

MP_integral_with_absolute_and_hypergeom_beta_(2).pdf