The zip command is used to apply a two-parameter procedure to two matched arrays (or lists) of arguments.

restart:
B:= (x,y)-> fsolve(y = (1-exp(-x*b))/(1-exp(-50*b)), b):
X:= #Some array of x values
Y:= #Some array of corresponding y values
zip(B, X, Y);

Bisection is a method of last resort. Why do you want to use it? You can quickly get any number of roots of f(p) with RootFinding:-NextZero. After making Preben's unapply correction, do this:

N:= 100:  #Number of roots desired.
Root:= Array(1..N):
Root[1]:= -RootFinding:-NextZero(p-> f(-p), 0):
for n to N-1 do
     Root[n+1]:= -RootFinding:-NextZero(p-> f(-p), -Root[n])
end do:

The minus signs are because it seems that you want the roots in the negative direction.

All that you need is

series(y(x+h) - y(x-h), h=0, 3);

NumberToList:= proc(x::realcons)
uses S= StringTools;
     parse~(S:-Explode(S:-Select(S:-IsDigit, sprintf("%a", x))))
end proc:

Change the last line of your last for loop from

theta[ii](tau):=rhs(%);

to

theta[ii](tau):= value(rhs(%));

The value command makes Maple perform the integrations.

Since the desired simplification is not generally true for all bases, you have to use the symbolic option to simplify:

simplify(%, symbolic);

Using Maple 17, when I use solve followed by evalf, I get

Try doing a restart, then try the solve and evalf again.

A real answer can be obtained with fsolve or RootFinding:-NextZero in addition to the method that you mentioned.

KernelDensityPlot is a command in the Statistics package that will output the same plot in the normal inline manner:

Statistics:-KernelDensityPlot(Cdata);

AB,C:= LinearAlgebra:-GenerateMatrix(
     [F[k] $ k= 0..3],
     [u[k,n] $ k= 0..3, u[k,n-1] $ k= 0..3],
     augmented= false
);
A:= AB[.., 1..4];  B:= AB[.., 5..8];

1. In most cases, you reduce the error by increasing the value of Digits, which has the default value 10.

Digits:= 15:
fsolve(eqs2);

2. To get an exact solution to a decimal problem (if it's possible), use convert(..., rational):

solve(convert(eqs2, rational));

The midpoint method is not implemented in dsolve(..., numeric, method= classical[...]), so you'll have to write your own. Here's my ad hoc version---tailored to this specific problem. As you can see, the error plot is a straight (in log-log mode) line with a slope of nearly exactly 2.

Download Midpoint_error.mw

You need to change cos(kz - `ωt`) to cos(k*z - omega*t). I ran the code, and everything works if you correct this expression in the three places that it occurs.

I am not sure what you mean by "the scale of plot". Do you want a logarithmic scale? Do you want to extend the range on the y-axes so that both plots have the same y-axis?

Your "final result" has three terms with dimensions of force and one term with dimensions of force/temperature. Was there an integral in your code that you expected would cancel the temperature, which is your variables alfa1 and alfa2? If so, I couldn't find it.

Kitonum wrote:

But you wrote that your method does not directly prove the uniqueness of the solutions found. And how much time require such proof?

Unique means one solution. So I would not use the word unique in this situation. We seek to prove that the two stated solutions are the complete set of solutions. When a solution is unique, my program will indeed prove it unique. (It is not absolutely guaranteed that it will always be able to complete such a proof; however, I've never seen a case with a unique solution where it failed to complete the proof.) So, to prove that the multiple solutions are the complete solution set, we add constraints so that we get proven unique solutions. If we add a complete (covers all logical possibilities) and mutally exclusive set of constraints (only one at a time, of course), and each addition produces (provably) 0 or 1 solutions, then we have all the solutions.

The additional time for my program to do the proof is about 0.02 seconds. Here it is:

Download City_LP.mw

Hmm. Somehow your default plot options got changed. Try this:

plot3d(x^2-y^2, x = -1 .. 1, y = -1 .. 1, style= patch);

If that shows the grid, then set the option as default by

plots:-setoptions3d(style= patch);

But I still wonder how your default got changed.