@Kitonum The unimodality is now proven (see above). There are at least two ways to make the fsolve work for large k. A simplistic way is to increase the value of Digits. Digits = 30 works for the k=50 case that you presented. A better way is to divide out a large positive factor, which, of course, doesn't change the zeros. It turns out that C(n,k) can be factored out of its own derivative with respect to k.

C:= (n,k,a,b)-> binomial(n,k)*a^(n-k)*b^k:
fsolve(simplify(D[2](C)(n,50,1,2)/C(n,50,1,2)), n= 50..infinity);

74.7487626691704

which agrees with your result.

@Axel Vogt Thanks. I have been thinking about this unimodality proof for a few days. To generalize, the generalized binomial coefficient function

C:= (n,k,a,b)-> binomial(n,k)*a^(n-k)*b^k

is unimodal as a function of k for 0 <= k <= n when a > 0 and b > 0. The derivative with respect to k is remarkably simple:

D[2](C)(n,k,a,b) = (Psi(n-k+1) - Psi(k+1) + ln(b/a))*C(n,k,a,b)

We can ignore the factor C(n,k,a,b) because it is strictly positive in the region of interest. We can ignore the shift by ln(a/b) because it does not depend on k. It remains prove that Psi(n-k+1) - Psi(k+1) is monotonic for 0 <= k <= n, which Axel did.

@oldstudent You can't use the plot command to plot equations. There are two options. The first is to subtract one side of the equation from the other and plot that, so that the x intercepts of the resulting plot are the roots of the equation:

plot(2*sin(2*x*Pi/180) - 1, x= 0..360);

The second is to plot separately the left side and right side of the equation, so that the x coordinates of the intersection points are the roots of the equation:

plot([2*sin(2*x*Pi/180), 1], x= 0..360);

You should get used to working in radians instead of degrees. All higher mathematics is done in radians, and nearly all mathematical software uses radians by default.

@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));