Your first version of eqd has a coefficient of omega0^2 on the y(t) term. The coefficient is omitted in the second version. If I include the coefficient in the second version, then I get the same phase portrait, the nearly concentric circles.

plot([B(t), A(t), t= t1..t2]);

where t1..t2 is some definite range of t-values.

For your first case, the derivative of f(x,y) with respect to t is, of course, 0. Perhaps you mean Diff(f(x,y,t), t, t)?

Your second case is a partial differential equation because the unknown function has more than one independent variable. Consequently, the command is pdsolve rather than dsolve.

pdsolve((Diff(f(x,y), (x,y), (x,y)))*x+(Diff(f(x,y), (x,y)))*y = 0, f(x,y));

In Kitonum's answer, take the line with the map command, and change the last character of that line from a colon (:) to a semicolon (;). Then you will see the output.

You could also do f(0.5), f(0.9), etc.

I wasn't chastising you for bad vocabulary. I truly didn't understand which alternative you meant. But Kitonum figured it out.

When you assigned sysdiff and fcns, you used curly braces ({ }). In your dsolve command you used them again. So the first argument is now a set containing a set rather than the set containing equations that is expected. So, changing the dsolve command to dsolve(sysdiff, fcns) will produce a solution. But like your other systems from today, the solution contains spurious constants of integration. So once again, make the command

dsolve(sysdiff, fcns, method= laplace);

This is very similar to an earlier problem you had from today where Preben suggested adding method= laplace to the dsolve arguments. That works here also, and it is the only thing that you need to add to get the plots. I don't know why it works. The reason that the result from the plain dsolve would not plot was that there were still integration constants left in the solution. (Should it be considered a bug?)

Remember, if you use dsolve(..., numeric) then you should use plots:-odeplot for the plotting (it is possible to use plot also with numeric, but it requires a different syntax).

It's essentially the same situation as with your last question. The difference is that this is solve rather than dsolve. You assign the result of the solving command to a variable. In this case, you used solution as the variable. Then my recommended technique is to access the contents of solution with the eval. To get the plot in this case, try

plot(eval(C2, solution), omega= 0..2, -3..3);

You wrote: How to get all solution?

I would think that with the huge number of similar Questions that you've asked (and have been Answered) over the past year that you'd be able to figure that out on your own by now.

I think that it is silly to ask for all solutions to this problem because there are a huge number, and most subsets of four points are a solution.

Nonetheless, all the solutions can be obtained with this code, which takes more than 30 minutes to run:

Sols:= table():
for P4 in combinat:-choose(L,4) do
     ct:= 0; #count non rt triangles for this group of 4 pts.
     for P3 in combinat:-choose(P4,3) do
          geom3d:-triangle(
               ABC, 
               [seq](geom3d:-point(A||k, P3[k][]), k= 1..3)
          );
          if geom3d:-IsRightTriangle(ABC) then  break  end if;
          ct:= ct+1
     end do;
     # 4 triangles can be made from 4 pts: C(4,3) = 4.
     if ct = 4 then  Sols[P4]:= [][]  end if 
end do:

Sols:= {indices}(Sols, nolist):
nops(Sols);

                             276392

A better technique is to use a Monte Carlo method to approximate the number of solutions. We use combinat:-randcomb to repeatedly select four points at random.

N:= 2^10:   #number of random trials
sol_ct:= 0: #number of solutions found
for k to N do
     P4:= combinat:-randcomb(L, 4);
     ct:= 0; #count non rt triangles for this group of 4 pts.
     for P3 in combinat:-choose(P4,3) do
          geom3d:-triangle(
               ABC,
               [seq](geom3d:-point(A||k, P3[k][]), k= 1..3)
          );
          if geom3d:-IsRightTriangle(ABC) then  break  end if;
          ct:= ct+1
     end do;
     # 4 triangles can be made from 4 pts: C(4,3) = 4.
     if ct = 4 then  sol_ct:= sol_ct+1  end if
end do:

trunc(evalf(sol_ct/N)*binomial(nops(L), 4));

                             278336

There are two ways: You can use dsolve(..., numeric) and plots:-odeplot, or you can use an analytic dsolve solution with a parametric plot.

Download phaseportrait.mw

I've experimented computationally with f(n) and discovered some very interesting properties. These would need to be proven of course.

Conjecture 1: If n > 2 and 2*n - 1 is prime, then f(n) is divisible by 2*n - 1.

I've verified this up n = 6940.

My second conjecture implies the first, and is much stronger. I've also verified it up to n = 6940.

Conjecture 2:If n > 2 and 2*n - 1 is prime, then 2*n - 1 is the largest prime factor of f(n).

So, even though f(6940) is a huge number with 4173 decimal digits, its largest prime factor is only 2*6940 - 1 = 13,879 (with only 5 digits). That's amazing!

Considering the huge size of f(n), my next conjecture "almost implies" the first two.

Conjecture 3: If n > 0, then f(n) has no prime factor greater than 2*n - 1.

Since your worksheet is long, I stopped reading after I found the first mistake. So, I don't know if this mistake is the ultimate cause of your issue. Correct the mistake, and if you still have your issue, then repost the corrected worksheet, and I'll have another look. (And if you don't still have the issue, please post a followup and say so!)

The mistake is that at one point you mispelled unapply as unpply. It's when you're defining I[Lm2] almost halfway through the worksheet.

You have three problems. The first is a trivial typo. Before the fprintf, you have a "*" which should be a ";". The second is that solve returns multiple solutions, some of which are not real. Your formatting is expecting a single real value. To get a single real value, use fsolve instead of solve (in both cases). The third problem is that with your current settings, each loop is going to execute 3 million times. Do you really want to spend the time required for that? If not, then you either need to reduce the 81920 or increase the value of dtt.

Download 1stMapleCalcI.mw

The desired expansion is not true for general A and B. (Consider A = B = -1, a = 1/2.) It is necessary to assume that they are positive. The sign of a doesn't matter.

Here are three ways:

expand((A*B)^a) assuming A>0, B>0;

expand((A*B)^a) assuming positive;

simplify((A*B)^a, symbolic);

[Edit: I completely redid this Answer, using all the original variable names, minus the underscores.]

There is a unique solution to your system. To reduce the visual space needed to present the solution, I removed the underscores from all of your names.

To be completely rigorous, we need to show that there is at least one assignment of numeric values to the coefficients that makes the denominator non-zero. This is of course trivial to do, and I leave it to you.

Download 16eqns.mw