@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

@reemeaaaah I have a habit of editing my answers a little too often. I just added some lines.
My knowledge of statistics is basically confined to the ability of spelling the word, in two languages, though.

@reemeaaaah I have a habit of editing my answers a little too often. I just added some lines.
My knowledge of statistics is basically confined to the ability of spelling the word, in two languages, though.

What effect do the extra arguments 100 and confidence=0.95 in the assignment
R := RandomTools:-Generate(
      distribution(Normal(-0.2e-1, 0.4e-1)),100,confidence=0.95, makeproc = true);
have?
It seems to me that they are just ignored. Try e.g. misspelling them or changing them in another way.
Before and after try
R();

The following syntax generates a list of 10 numbers
R := RandomTools:-Generate(list(distribution(Normal(-0.2e-1, 0.4e-1)),10));

and this
R := RandomTools:-Generate(list(distribution(Normal(-0.2e-1, 0.4e-1)),10),makeproc=true);
makes a procedure which at each invocation does the same thing:
R();

I was writing in frustration after having answered and afterwards edited the answer to:

http://www.mapleprimes.com/questions/143595-E-With-One-Random-Variable

At least my frustration got vented.

I was writing in frustration after having answered and afterwards edited the answer to:

http://www.mapleprimes.com/questions/143595-E-With-One-Random-Variable

At least my frustration got vented.

@reemeaaaah The MaplePrimes editor very often leaves out everything on a line after <.
That is quite a nuissance.
I corrected it. I hope it is OK now.

@reemeaaaah The MaplePrimes editor very often leaves out everything on a line after <.
That is quite a nuissance.
I corrected it. I hope it is OK now.

Well, I did look at it.

Did you try

evalf[20](eval(yg,p=0.3456));

evalf[20](eval(wronski,p=0.3456));

There is a huge discrepancy. My guess is that your wronski is wrong. I didn't try to follow the argument.

At closer inspection I now see that the angular frequencies appearing in sin and cos are v and 2*v, where
v = 0.5995250000.

I start as before with getting rid of that extra beta in beta(.787*beta). It should not have been there in the first place. I wonder how it got in there?

My guess is that the floats appearing hide a simple symbolic relationship beyond what I have already mentioned.
So we shall try to find that.

restart;
Digits:=10:
A := Int((1.334389725*sin(1.199050000*beta)/beta+(1.669308692*(cos(1.199050000*beta)-1.667987156*sin(1.199050000*beta)/beta+2.782181154*sin(.5995250000*beta)^2/beta^2))/beta^2)^2/((61.09*(240.38915*tanh(.787*beta)*beta(.787*beta)+1))/(61.09+3.935*tanh(.787*beta)*beta(.787*beta))+1.7314*coth(.787*beta)*beta(.787*beta)), beta = 0 .. infinity);

indets(A,specfunc(anything,beta));
subs(beta(.787*beta) = beta*(.787*beta) , A);
B := normal(%);
f := IntegrationTools:-GetIntegrand(B);
indets(f,specfunc(anything,{sin,cos}));
subs(1.199050000*beta=2*v*beta,.5995250000*beta=v*beta,f);
op([2,1],%);
Cs:=coeffs(%,[beta,cos(2*beta*v),sin(2*beta*v),sin(beta*v)],'T');
T;
w:=add(a[i]*T[i],i=1..nops([T]));
param:= convert(zip(`=`,[seq(a[i],i=1..nops([T]))],[Cs]),set);
#For the integral to be convergent on the interval 0..1 it is necessary and sufficient that the series expansion:
series(w,beta=0,5);
#contains no powers of beta under 4:
eval(%,param union {v=.5995250000});

It so appears. But rounding as seen in my earlier response can very easily spoil the symbolic relationship:

a[2]+2*a[3]*v+a[4]*v^2 = 0

My advice is to keep symbolic input until you are forced to resort to numerics.

At closer inspection I now see that the angular frequencies appearing in sin and cos are v and 2*v, where
v = 0.5995250000.

I start as before with getting rid of that extra beta in beta(.787*beta). It should not have been there in the first place. I wonder how it got in there?

My guess is that the floats appearing hide a simple symbolic relationship beyond what I have already mentioned.
So we shall try to find that.

restart;
Digits:=10:
A := Int((1.334389725*sin(1.199050000*beta)/beta+(1.669308692*(cos(1.199050000*beta)-1.667987156*sin(1.199050000*beta)/beta+2.782181154*sin(.5995250000*beta)^2/beta^2))/beta^2)^2/((61.09*(240.38915*tanh(.787*beta)*beta(.787*beta)+1))/(61.09+3.935*tanh(.787*beta)*beta(.787*beta))+1.7314*coth(.787*beta)*beta(.787*beta)), beta = 0 .. infinity);

indets(A,specfunc(anything,beta));
subs(beta(.787*beta) = beta*(.787*beta) , A);
B := normal(%);
f := IntegrationTools:-GetIntegrand(B);
indets(f,specfunc(anything,{sin,cos}));
subs(1.199050000*beta=2*v*beta,.5995250000*beta=v*beta,f);
op([2,1],%);
Cs:=coeffs(%,[beta,cos(2*beta*v),sin(2*beta*v),sin(beta*v)],'T');
T;
w:=add(a[i]*T[i],i=1..nops([T]));
param:= convert(zip(`=`,[seq(a[i],i=1..nops([T]))],[Cs]),set);
#For the integral to be convergent on the interval 0..1 it is necessary and sufficient that the series expansion:
series(w,beta=0,5);
#contains no powers of beta under 4:
eval(%,param union {v=.5995250000});

It so appears. But rounding as seen in my earlier response can very easily spoil the symbolic relationship:

a[2]+2*a[3]*v+a[4]*v^2 = 0

My advice is to keep symbolic input until you are forced to resort to numerics.

@mtj4u 

This time the function beta(.787*beta) appears. Does it mean (1) .787*beta, (2) beta, or (3) beta*(.787*beta)?
In any case f will have a singularity at zero sufficient to make the integral from 0 to 1 divergent.

restart;
Digits:=100:
A := Int((1.334389725*sin(1.199050000*beta)/beta+(1.669308692*(cos(1.199050000*beta)-1.667987156*sin(1.199050000*beta)/beta+2.782181154*sin(.5995250000*beta)^2/beta^2))/beta^2)^2/((61.09*(240.38915*tanh(.787*beta)*beta(.787*beta)+1))/(61.09+3.935*tanh(.787*beta)*beta(.787*beta))+1.7314*coth(.787*beta)*beta(.787*beta)), beta = 0 .. infinity);

indets(A,specfunc(anything,beta));
subs(beta(.787*beta) = beta*(.787*beta) , A); #First interpretation
B := normal(%);
f := IntegrationTools:-GetIntegrand(B);
series(f, beta = 0, 9);

@mtj4u 

This time the function beta(.787*beta) appears. Does it mean (1) .787*beta, (2) beta, or (3) beta*(.787*beta)?
In any case f will have a singularity at zero sufficient to make the integral from 0 to 1 divergent.

restart;
Digits:=100:
A := Int((1.334389725*sin(1.199050000*beta)/beta+(1.669308692*(cos(1.199050000*beta)-1.667987156*sin(1.199050000*beta)/beta+2.782181154*sin(.5995250000*beta)^2/beta^2))/beta^2)^2/((61.09*(240.38915*tanh(.787*beta)*beta(.787*beta)+1))/(61.09+3.935*tanh(.787*beta)*beta(.787*beta))+1.7314*coth(.787*beta)*beta(.787*beta)), beta = 0 .. infinity);

indets(A,specfunc(anything,beta));
subs(beta(.787*beta) = beta*(.787*beta) , A); #First interpretation
B := normal(%);
f := IntegrationTools:-GetIntegrand(B);
series(f, beta = 0, 9);

This time you somehow introduced another beta function: beta(.787*beta).
So the question is: what was the intention? Is it supposed to mean just  .787*beta, or is it supposed to mean beta, or beta*(.787*beta)?

In all 3 cases f has a serious singularity at zero making the integral over beta = 0.. 1 divergent.

restart;
Digits:=100: #overdoing it some maybe.
A := Int((1.334389725*sin(1.199050000*beta)/beta+(1.669308692*(cos(1.199050000*beta)-1.667987156*sin(1.199050000*beta)/beta+2.782181154*sin(.5995250000*beta)^2/beta^2))/beta^2)^2/((61.09*(240.38915*tanh(.787*beta)*beta(.787*beta)+1))/(61.09+3.935*tanh(.787*beta)*beta(.787*beta))+1.7314*coth(.787*beta)*beta(.787*beta)), beta = 0 .. infinity);

indets(A,specfunc(anything,beta));
subs(beta(.787*beta) = .787*beta, A); #Taking the first interpretation mentioned
B := normal(%);
f := IntegrationTools:-GetIntegrand(B);

asympt(f, beta, 9); #May not work in Maple 14 (but irrelevant)
series(f, beta = 0, 9);

My preference would be

evalindets(sol,specfunc(anything,{sin,cos}),s->op(0,s)(omega*t));

But your way could be made similarly:

evalindets(sol,specfunc(anything,{sin,cos}),x->applyrule(eq,x));