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.

Actually, the representation of the table as matrices fits across my screen.

restart:
interface(rtablesize= 27):
M:= [seq(seq(seq(< < a,0 > | < b,c > >, a= 0..2), b= 0..2), c= 0..2)]:
T:= Matrix(27,27, (i,j)-> M[i].M[j] mod 3);

The piecewise command is easier to use than a chain of if-then-else statements.

Grades:= proc(x::realcons)
local a,b,c,d,f;
     piecewise(
          x < 0, undefined, x < 59.5, f, x < 69.5, d, x < 79.5, c,
          x < 89.5, b, x <= 100, a, undefined
     )
end proc:

If you still want to use if statements, then use elif, which is short for else if. That way, you can turn a chain of nested if statements into a single statement:

Grades:= proc(x::realcons)
local a,b,c,d,f;
     if x < 0 then undefined
     elif x < 59.5 then f
     elif x < 69.5 then d
     elif x < 79.5 then c
     elif x < 89.5 then b
     elif x <= 100 then a
     else undefined
     end if
end proc:

You need a for loop, an if-then statement, the isprime command, and a counter variable.

CountPrimes:= proc(n::posint)
local j, count:= 0;
     for j from 2 to n do
          if isprime(j) then
               count:= count+1
          end if
     end do;
     count
end proc:

CountPrimes(15);
                               6

Go to the Maple Applications Center, and download and install the OrthogonalExpansions package. Then

OrthogonalExpansions:-FourierSeries(evalc(x*exp(I*x)), x= -Pi..Pi, infinity);

The evalc is necessary to get the correct coefficients at i=1. Don't confuse the summation index i with the imaginary I.

The last constraint, constraint 21, "Eilish is not next to the person with straight hair", needs to be specified as Rel(Separated, Eilish, straight, PN, [1]), which says that they are separated by at least 1 when measured wrt position. You merely have Eilish <> straight. Since Separated is an export of the dynamic module, this constraint needs to be specified after the with of the module.

When I do this, the program reports multiple solutions. The second solution that it gives is the same as the solution on the website from which you got the problem. From my reading, the first solution that the program gives also satisfies every constraint. Please verify this.

Here's how to get the answers with fsolve. This can be modified to instead use the multivariate Newton procedure that Markiyan directed you to.

restart:
eqns:= {
     t*((p-.764*z-2.194768)^2-1.170308549+.529948*(z-.382)^2)-(1-t)*p,
     t*((p+.382*z+.661624*y-1.907568)^2-1.018097144+.529984*(-.866*y-.5*z-.382)^2)+(1-t)*y,
     t*((p+.382*z-.661624*y-2.470348)^2-1.31636154+.529984*(.866*y-.5*z-.382)^2)+(1-t)*z
}:

(a,b):= (0,1):  h:= 0.0001:
N:= round((b-a)/h)+1;

Sols:= Array(1..N, 1..4):
inits:= {y=1,z=1,p=1}:
tt:= a:
for n to N do
     Sol:= fsolve(eval(eqns, t= tt), inits);
     Sol:= eval([y,z,p], Sol);
     Sols[n,..]:= < tt, Sol[] >;
     tt:= tt+h;
     inits:= {([y,z,p] =~ Sol)[]}
end do:

The solutions can be plotted as a space curve as a function of t like this:

plots:-spacecurve(convert(Sols[.., 2..4], listlist), labels= [y,z,p]);

Use an animation, so that one of a, c, n becomes the time parameter. Here I use a as the time parameter.

restart:
f:= (a,c,n)->
     64*a^2*n-16*a^4*n-256*c^2*n+32*c^4*n-96*a^2*c^2+ 8*a^2*c^4+24*a^4*c^2
     -64*a^2*n^2+256*a^2-124*a^4+ 15*a^6+256*n^2+64*a^2*c^2*n
:
plots:-animate(plot3d, [f(a,c,n), c= 0..2, n= 0..4], a= 0..2);

The above ideas can be made into a custom convert procedure:

`convert/decimalfraction`:= (x::float)-> op(1,x)/``(10^(-op(2,x))):

convert(.25, decimalfraction);

x^(1/3) is complex valued for negative x because Maple uses the principal branch with `^`. The real-valued alternative is surd.

plot(surd(x,3), thickness= 5);

You can generate the procedure by calling LinearAlgebra:-Eigenvectors on a symbolic 2x2 matrix A = < < a, b > | < c, d > > and then reverse engineering the results. Doing that, I got the following

Eigen:= proc(A::'Matrix'(2,2))
local
     a:= A[1,1], b:= A[2,1], c:= A[1,2], d:= A[2,2],
     s:= sqrt((a-d)^2+4*b*c)/2,
     e1:= (a+d)/2 + s, e2:= (a+d)/2 - s,
     ev1:= < c, (d-a)/2+s >, ev2:= < c, (d-a)/2-s >
;
     if s=0 then [] else [e1, e2, ev1, ev2] end if
end proc:

This procedure needs a little modification. There are cases where s=0 but there are still two linearly independent eigenvectors. For example, any nonzero vector is an eigenvector of the zero matrix.