@matja There must have been something else interfering with the assume command since the following works:

restart;
assume(a>0);
limit(exp(-a*x), x = infinity);

Certainly S1+S2 is the sum of the two expressions, so I don't know what your problem is. Do you want to use simplify, expand, or combine on the sum or what?

An animation:

plots:-animate(plot, [S, theta = 0 .. 10], t = 0 .. 1, frames = 100);

The sum in S1 is the first few terms of

Sum(cos((2*k+1)*theta)/(2*k+1)*(-1)^k,k=0..infinity);

which is the Fourier series of

f:=piecewise(abs(theta)<Pi/2,Pi/4,-Pi/4);

restart;
eq1:=eval(x*tan(x)-I*i*Pi/(1+I*i*Pi), i = 10);
eq2:=eval((1+I*i*Pi)*x*sin(x)-I*i*Pi*cos(x), i = 10);
res:=RootFinding:-Analytic(eq2, x = -20-I .. 20+I);
nops([res]);
#Check:
evalf(eval~(eq1,x=~[res]));
simplify(fnormal(%,8));

restart;
eq1:=eval(x*tan(x)-I*i*Pi/(1+I*i*Pi), i = 10);
eq2:=eval((1+I*i*Pi)*x*sin(x)-I*i*Pi*cos(x), i = 10);
res:=RootFinding:-Analytic(eq2, x = -20-I .. 20+I);
nops([res]);
#Check:
evalf(eval~(eq1,x=~[res]));
simplify(fnormal(%,8));

@Markiyan Hirnyk By doing successively

RootFinding:-Analytic(eval(x*tan(x)-I*i*Pi/(1+I*i*Pi), i = 10), x = -20-I .. 20+I);
RootFinding:-Analytic(eval(x*tan(x)-I*i*Pi/(1+I*i*Pi), i = 10), x = -20-I .. 20+I,continue);

you find 6 roots.

@Markiyan Hirnyk By doing successively

RootFinding:-Analytic(eval(x*tan(x)-I*i*Pi/(1+I*i*Pi), i = 10), x = -20-I .. 20+I);
RootFinding:-Analytic(eval(x*tan(x)-I*i*Pi/(1+I*i*Pi), i = 10), x = -20-I .. 20+I,continue);

you find 6 roots.

You are more likely to get an answer if you give us the text instead of an image.

Which parameter values will make it run? 
Maybe an uploaded worksheet will help.

Which parameter values will make it run? 
Maybe an uploaded worksheet will help.

Here is a version which I believe is sound. It adds and subtracts series expansions of sin and cos to make the individual integrals convergent. It is essential that the sum of the terms in L1 be 0. That follows from knowing that the original integrand has no singularity at zero.

The code should work for even as well as odd values of N.
It is still slow for all but one digit values of N.

restart;
with(IntegrationTools):
A:=n->Int(mul(sin(x/2^k)*2^k/x, k = 0 .. n), x = 0 .. infinity);
#N:=9:
N:=4:
Change(A(N),x=2^N*t):
B:=combine(%):
f1:=GetIntegrand(B):
series(numer(f1),t=0,N+1); #terms of order up to and including N are zero
L:=convert(expand(f1,sin,cos),list):
sn:=unapply(convert(series(sin(x),x=0,N+1),polynom),x):
cs:=unapply(convert(series(cos(x),x=0,N), polynom),x):
L1:=eval(L,{sin = sn,cos=cs}):
simplify(convert(L1,`+`)); #Will return 0.
map(int,L-L1,t=0..infinity):
convert(%,`+`);

Here is a version which I believe is sound. It adds and subtracts series expansions of sin and cos to make the individual integrals convergent. It is essential that the sum of the terms in L1 be 0. That follows from knowing that the original integrand has no singularity at zero.

The code should work for even as well as odd values of N.
It is still slow for all but one digit values of N.

restart;
with(IntegrationTools):
A:=n->Int(mul(sin(x/2^k)*2^k/x, k = 0 .. n), x = 0 .. infinity);
#N:=9:
N:=4:
Change(A(N),x=2^N*t):
B:=combine(%):
f1:=GetIntegrand(B):
series(numer(f1),t=0,N+1); #terms of order up to and including N are zero
L:=convert(expand(f1,sin,cos),list):
sn:=unapply(convert(series(sin(x),x=0,N+1),polynom),x):
cs:=unapply(convert(series(cos(x),x=0,N), polynom),x):
L1:=eval(L,{sin = sn,cos=cs}):
simplify(convert(L1,`+`)); #Will return 0.
map(int,L-L1,t=0..infinity):
convert(%,`+`);

@Markiyan Hirnyk I overlooked the fact that the sine integrals in the first part of my answer are also divergent.
I have inserted notes in my answer.
I prefer correct reasoning with a wrong result over wrong reasoning with a correct result.

@Markiyan Hirnyk I overlooked the fact that the sine integrals in the first part of my answer are also divergent.
I have inserted notes in my answer.
I prefer correct reasoning with a wrong result over wrong reasoning with a correct result.

Have you managed to make this work in a simpler case?
For instance with r0 =1  and rf = 10 and number of partitions = 2?

@reemeaaaah In order to have total control over the function calls I thought it safer to use a bare bone version of Euler's method.
I made two versions. Euler is the usual and requires input of the function f in y' = f(t,y).
Euler2 requires input of a function f of 3 variables y' = f(t,y,r), where the third argument is meant for the random input.
Your particular normal distribution is hardcoded into Euler2.
For what it is worth, here it is:

StochasticEuler.mw