@zenterix Here is what I got when running your file on my Linux machine. 

In the first attempt Maple crashed after running for some half an hour or so.  I montitoried the system's status during the execution and observed that my swap partition, which is a measly 0.8 GB, got filled and apparently that was the cause of the crash.

Then I added a largish (100 GB) swap file and ran the program again.  I didn't time it, but I am guessing that it took about an hour or so to produce the sphere with an arrow.  It turned out that 100 GB of swap space was much more than necessary.  Only a small fraction of it was used.  As I played the animation, however, the memory usage went up until it nearly ate up all of the 64 GB RAM of my computer's memory and the machine started swapping the memory to disk.   At that point I killed the program because the machine had become almost unresponsive.

With your 256GB RAM you may have better luck. 

The lesson learned here is you need a descent sized swap space to run your program.  Type "free" on your Linux command-line to see the size of your swap area.  (The number shown is in kilobytes).  It looks like you will need something like 20 GB of swap space to run your program.

Normally a swap partition is created during the installation of Linux.  If your swap partition is too small, you may temporarily add a swap file to enlarge it.  It's not difficult.  Let me know if you need help with doing that.  Perhaps your Mac machine has a large swap partition and that's why things work smoothly there.

@Carl Love Your post came just in the right time. I was looking for a way for doing just that.  Vote up!

In an attempt to understand what you are asking, I searched for "Bart Simpson, epicycles" on Google and found a YouTube Video which I think is related to your question.  Let me summarize.

We have a line drawing in the complex plane.  We parametrize the drawing as F(t), where t goes from 0 to 2*Pi.  Then we expand F(t) into the complex Fourier series, and interpret the partial sums of the series as depicting a planet, a moon, a moon of that moon, and so on.  Thus, the Fourier series approximation of F(t) may be viewed as the orbit of the last moon in the sequence of the moons.

In the attached Maple worksheet I show how to animate the moons as they trace the curve corresponding to F(t), which, for the sake of illustration I have taken to be a square.  The animation below corresponds to picking 7 terms of the Fourier series, and consequently the target square (drawn in green) is rendered with rounded corners (rendered in red).

Download worksheet: mw.mw

You may replace the square with any other drawing provided that you can obtain a parametric representation for it.  I understand that you are looking for a method of extracting the data for a given curve from an image, and then parametrizing it.  That's a rather complex operation.  I haven't given it much thought.  The video quoted above suggests using a traveling salesman algorithm to connect points extracted from a graph.

Added Dec 5:

Here I have replaced the square with an F-shape:

@ecterrab I updated to the latest version of the Physics package and the trouble went away.  I have added a postscript to my earlier post indicating that pdsolve() was not at fault.

@ecterrab Hi Edgardo, I have narrowed down the problem to the subscripted arbitraryfunctions business.

Specifically, in my mapleinit file I have:
`pdsolve/arbitraryfunctions` := "subscripted":

If I remove that line, then 
pdsolve({pde1, pde2, ic1, ic2});
works just fine.  Otherwise I get the error that I had noted.

The attached worksheet was executed with the troublesome line commented out in mapleinit and it shows that pdsolve works in that case.  Further down in the worksheet I insert
`pdsolve/arbitraryfunctions` := "subscripted":
after which we see that pdsolve fails.

Download mw2.mw

@astroverted My original worksheet works in Maple 2021 and 2022.  It fails in earlier versions.  Here is an updated worksheet that works in Maple 2017 onward.

In this worksheet I have set the grid sizes m and n to small values to speed up the calculations.  Increase the values for greater accuracy (but then you will have to be patient while you wait for the results).

Download strange-PDE-ver2.mw 

@imparter You're welcome.  I have reported this bug to Maplesoft and I hope that it will be corrected in the next release.

@astroverted 

dharr's assertion regarding the boundary conditions is correct.  The paper that you have cited provides four initial/boundary conditions in equations (15) through (18).  But (16) and (18) are not arbitrary; they are derived from (15) and (17) and the PDEs as I will explain in the Answer that I am going to post next.  The paper assumes that's obvious and does not state that fact explicitly.

@romanrieme As acer suggested, you need display(F2(t),F3(t)):

animate(display, [`if`(t < 0, F2(t), display(F1(t), F3(t))), scaling = constrained], t = -1 .. 1, frames = 140);

Alternatively, you may do it as I suggested:

display([animate(F2, [t], t = -1 .. 0), animate([F1, F3], [t], t = 0 .. 1)], insequence)

@mmcdara I use "extrapolation" for that purpose, as in "φ(x,y) extrapolates the boundary data into the rectangle".  I have seen people refer to that as "Dirichlet lift" or just "lift", which is closer to the French "relèvement".

You can expect better feedback in this forum if you upload your worksheet rather than what you have posted.  But just looking at the first few lines of your code, I see that you have strange stuff such as

ibvc:=u(x, 0)=e^(K*x)/(1+exp^(0.5*K*x))^2,v(x, 0)=1/(1+exp^(0.5*K*x));

What is e^(K*x)?   What is exp^(0.5*K*x))?  Think!

@mmcdara Yes, with u(x,0)=0 it's possible to solve the problem in Maple by Fourier series.

Here is the outline of the idea.

We reflect the domain about the line y=0 and thus obtain a rectangle.  We solve Laplace's equation on the rectangle by taking for boundary values on the new edges the negatives of the values of the opposite pages.  Then the solution will be zero automatically along the line y=0 due to symmetry.

A detail to be taken care of is that the rectangle's edges are not parallel to the coordinate axes.  To apply the Fourier series, we rotate the rectangle by 45 degrees to bring its edges into alignment with the coordinate axes.  The PDE does not change since the Laplacian operator is rotationally invariant.

To avoid the rotation step, it would be much easier to rotate the original domain.  That is, instead of the triangle shown in the OP's drawing, we take a triangle formed by the coordinate axes and the line x+y=1 (or ax+by=c for some constants a, b, c).  Then the reflection of the triangle across the hypotenuse will produce a rectangle in which we may apply the standard Fourier series in the straightforward way.

Even then, there is an extra snag that needs to be taken care of.  The Fourier series method requires homogeneous boundary conditions, that is, u=0, along the four edges of the rectangle.  But our boundary data is not zero along the edges.  To take care of that, we introduce a function φ(x,y)  defined over the rectangle so that φ matches the given boundary conditions along the rectangle's edges.  If we let w(x,y) = u(x,y) - φ(x,y),  then w will be zero over the rectangle's edges, and therefore Fourier series may be applied to calculate it.  We note, however, since u satisfies the Laplace's equation, w satisfies the Poisson equation w_xx + w_yy = −φ_xx −φ_yy. which may be solved by Fourier series just as easily.

Finally, there remains the question of how one constructs the function φ(x,y).  That is the subject of a topic that I posted here a few years ago.

There is some work involved in implementing the steps outlined above.  I haven't bothered to do the work considering that this approach applies to the very limited case when the domain is a right triangle, the boundary condition on the triangle's hypotenuse is zero, and also the fact that the OP has specified no such boundary condition.  In general such problems are best solved through the method of finite elements which is not available in Maple.  

@Oliveira Sure, you get a Fourier series representation of the solution.  That works on a rectangle but your original question was about solving over a triangle.

@bstuan Let f(x) be your integrand and let m be the minimum value of (x+2)/(2x^3+1) on the interval (1,2). Since f(x) > m/ln(x), the integral of f(x) is more than the integral of m/ln(x) on (1,2). But the integral of m/ln(x) on (1,2) is infinity, and therefore so is the integral of f(x).

Hint: Examine the convergence of the integral of 1/ln(x) over the interval (1,2).  That should be good enough to reach a conclusion about your original problem.