@janhardo It does appear that `mod` is defined in the .m file. To see it do

showstat(`mod`);

Note the backquotes, which are required since mod is a keyword.

There may be procedures inside modules in your file. The anames process won't find those. So do this also:

before:= {anames(`module`)}:
read "MyFile.m":
newmods:= {anames(`module`)} minus before;

Note once again the backquotes, which are required because module is a keyword.

If you do have modules in the .m file, let me know. Listing the procedures in modules is a little bit more complicated.

I think we need some type of general form for the final expression, such as a^i + b^j with i, j integer.

@janhardo Once procedures are stored in .m format, you cannot see them exactly as they were originally written. However, you can see something fairly close to the original using showstat. This still requires knowing the names of the procedures. So, if you know that there is a procedure named MyProc, do

showstat(MyProc);

@Markiyan Hirnyk I do use it because I have to know that the point returned by fsolve is the unique point of zero derivative.

@Markiyan Hirnyk I mean that it is unimodal when viewed as a function on the real numbers, as in this plot

plot(binomial(5,k), k= 0..5);

The more general question is How do you convert a "random variable" (as that term is used by the Statistics package) into a "stochastic process" (as that term is used by the Finance package).

I only see ODEs, not a PDE.

Anyway, just include the equation for s(t) with the others, provide a value for r, provide an initial condition for s(t), and run dsolve again.

@toandhsp See the Reply by Markiyan. The equation has two real solutions. Using RealDomain, the equation has one real solution. In my opinion, RealDomain is correct to reject the solution -1-sqrt(2). If RealDomain were simply a filter to discard non-real solutions, it would be rather trivial.

@dcanright Matrices are not just intended for algebra. This is just a little bug that prevents you from using a one-dimensional initializer. Using a two-dimensional initializer of strings works fine:

Matrix(2,2, [["A", "B"], ["C", "D"]]);

@Markiyan Hirnyk In the 3d plot that you show, one function is linear and the other is piecewise linear. To me, those are not believable representations of the functions.

@radzys Unfortunately, using your suggested fifth boundary condition (D@@3)(phi)(L) = (D@@2)(phi)(L) leads to the error message "Matrix is singular", probably because it leaves the possibility of an identically zero solution.

@Markiyan Hirnyk 

The collect(combine(expand(...))) does nothing is this case. Indeed, they leave p unchanged. Only the coeff(p, x^n)  (rather than the OP's coeff(p, x, n)) is needed. The following also work:

coeff(p, x^(n+2));

coeff(coeff(p, x^(n+2)), y^(n+2));

Your system contains parameters omega and k_r for which no numerical values are given. Do you have values for these?

Runge-Kutta or rkf45 can only be used for initial value problems (IVPs). This is a boundary value problem (BVP). Better to just let dsolve choose the method.

You have numerous syntax errors which I'll show you how to correct once you supply the values of the parameters.

@wo0olf 

Use the output of dsolve as a function:

sol1(.0001);

@hacker9130368 Going from intensity=9000000 to intensity=10000000, I do get an increase in b and a decrease in c.

With intensity=9000000:

With intensity=10000000:

If you got no difference, then you probably need to set Digits to a higher value. I used Digits:= 15.