I believe the solution below to be more computationally stable. It uses fsolve to avoid the issue of multiple complex solutions returned by solve. For the derivative of the implicit function, I used the standard formula for that derivative in terms of the partial derivatives of the implict equation: dy/dx = - (dF/dx) / (dF/dy) where F(y,x) = 0.

Download pH.mw

Here's a substantially different approach to parallelizing the code. This one does show a benefit to parallelizing, provided that you can accept the solution divided into several chunks. Being divided into chunks might be beneficial for subsequent parallelizing.

The trick in this technique is to divide the side vectors into equal-sized (or nearly equal-sized) chunks such that there is one chunk per processor. Then we use Threads:-Seq with the [tasksize= 1] option.

Download ParallelModular2.mw

I got a partial analytic solution and a numeric solution. There is something that I don't inderstand about the analytic solution. My question is in the worksheet. Hopefully someone else will answer this.

Download numPDEsol.mw

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.