@MarcinM The documentation says that plots:-spacecurve takes a list or an rtable as an argument. A Vector is a special case of the rtable, so we are okay.

The documentation of plot3d says that the input should be a list.  The use of a Vector as input does not seem to be documented but that has worked for as long back as I can remember.

@MarcinM Regarding "I do not understand the philosophy of plotting vector functions in Maple", you need to go back to Stewart's chapter 13 where space curves are introduced as vector-valued functions like this:
    r(t) = < f(t), g(t), h(t) >
where t is a parameter.

The Maple command spacecurve( < f(t), g(t), h(t) > , t=a..b) plots that parametric curve.

But the parameter need not be named "t".  Any symbol may be used for the parameter.  In the answer that I posted, the parameter is named "x", as in:
  C1 := < x, sqrt(1-x^2), 0>;

As to plotting surfaces in Maple, we do

plot3d([f(x,y), g(x,y), h(x,y)], x=a..b, y=c..d);

The parametric surface here is described as a list [f, g, h], but a vector description, as in <f,g,h> works just as well and produces the same result.  I used the vector representation in my answer.

@simplevn1967 You referred to "Maple 17".  Did you perhaps mean "Maple 2017"?  The two are quite different.   "Maple 17" was released in 2013 while  "Maple 2017" was released in 2017.

I ran the worksheet that you posted in my Maple 2017 and it worked just fine.  I don't have "Maple 17" on my machine so I cannot test it.  Perhaps "Maple 17" is too old for that code.

@simplevn1967 I cannot diagnose the problem without seeing your worksheet.

As you reply to this message, you will see a large green upward pointing arrow in the toolbar.  Click on it to attach your worksheet to your reply.

@mmcdara Your calculation imposes the initial condition y(0)=0 but the result does not satisfy y(0)=0.  In what way is that a solution to the initial value problem?

@Carl Love That's cool, but perhaps too clever for my taste.  I prefer more verbose code that one can verify at a glance rather a overly terse code that requires effort to decipher.   Vote up, nonetheless.

@mmcdara I had no intention of criticizing your suggestion.  My remark was directed to the OP, who in response to my question of whether he knew how to write the equation of a tangent plane to a sphere, he responded with an appeal to Maple's geom3d package, while the simple answer is just (X − P) . (P − C) = 0, where C is the sphere's center, P is the point of tangency, and X is an arbitrary point on the tangent plane.

@rlopez Yes, I recall having seen finite difference schemes applied to non-rectangular domains. I don't think that's worth the effort.

Writing a finite element scheme using 1st degree elements will take just a few dozens lines of Maple code and it will work for arbitrary-shaped domains. That's provided that we have a general proc for triangulating the domain. Maple's GraphTheory:-DelaunayTriangulation comes close to fulfilling that task but is not quite adequate.  For instance, it is not suitable for triangulating nonconvex domains.  Even on convex domains Delaunay triangulation produces poor results for finite element applications.

Jonathan Shewchuck's Triangle library (written in C) produces perfect triangulations of arbitrary two-dimensional domains bounded by polygons and possibly containing one or more polygonal holes.  I have made extensive use of that code in my programs (also in C).  I suspect that it would be possible to compile Shewchuck's C code into a Maple library but I haven't looked into that.  Once that's done, a slew of possibilities will open for implementing finite element code in Maple.  Perhaps someone at Maplesoft would be interested in looking into that.

It's not clear how much help you will need for this.  Do you know how to write the equation of the tangent plane to a sphere at a given point?  If you explain how you would do that, then I will be happy to show how to do that for a list of points.

@rlopez That's a nice video.  In it you explain a variety of ways, and their pros and cons, of solving elliptic and parabolic PDEs.

The finite differences method that lies at the foundation of all those methods, limits the spatial domain to a rectangle (of any dimension).  It will be difficult, if not impossible, to extend those methods to more general domains. (I am thinking of a square domain from which one or more holes have been punched out).  That's not a fault of your presentation; that's an inherent limitation of the finite difference methods.

At one point you say: "I have a feeling that this is why Maple does not have an elliptic numeric solver cause it is hard to have it robust enough to handle a variety of problems."

I agree.  Implementing a solver that works on arbitrary-shaped domains is a monumental task. There are specialized companies—Comsol and LS-Dyna come to mind—devoted to that task.  Their underlying method is finite elements rather than finite differences.  Perhaps one day Maple will branch out in that direction but I don't envision that happening anytime soon.

@sand15 If the predator starts at an arbitrary position, of course there is no reason for it to capture the target exactly at the origin.

The point of the puzzle is to find out where the predator is to start so that the capture occurs at the origin.  The answer is:  if the predator starts anywhere on an ellipse (that is to be determined), then the capture does occur at the origin.

@sand15 Yes, you are quite correct.  I had not explored the possibility of such small values of E1 and E2. Thanks for pointing that out.

@raj2018 I am getting zero for f1 evaluated at phi=0, as you expect.  Have a look.

Download N-eval-2.mw

You say that "the eval command does not gives the correct value of its first derivative" but your worksheet has no derivative in it.  Perhaps you meant something else?

@Muhammad Usman Floating point approximations can affect eigenvalue calculations severely.  I wouldn't recommend changing exact values to floating point approximations merely to speed up the calculations.  Running your worksheet without the evalfs with M=50 takes less than 8 minutes on my (very old) computer.  About half of that time is spent on calculating the matrix B.  The other half is spent on calculating the determinant.  Solving  the determinant for the eigenvalues takes a fraction of  a second; see the attached worksheet.

8 minutes is not quick, but it's worth the wait for having the assurance that the eigenvalues are computed correctly.

Worksheet: eigenvalues-50.mw