This looks like a homework problem.  Here I will provide most of the solution but in the future please show some evidence of what you can do on your own before appealing for help.

Download mw.mw

Download mw.mw

This worksheet shows how to plot those curves.  These do not agree with the picture that you have shown.  Perhaps the picture corresponds to a different choice of parameters.

The circle marks the solution of the system {eq1, eq2}.

Download mw.mw

I think the proper way to doing this is to work out the 2D to 2D transformation that deforms your equirectangular projection (I will refer to is as "the picture") onto the final result, without going through the intermediate 3D representation.

The following worksheet can serve as a starting point.  Let 0 < s < 2*Pi and -Pi/2 < t < Pi/.2 be the horizontal and vertical coordinates of the picture.  The mapping
(s, t) -> [cos(t)*cos(s), cos(t)*sin(s), sin(t)]
maps the picture into the sphere in 3D.  This is shown in plot p1.

The mapping
(s, t) -> [sqrt(2*(1+sin(t)))*cos(s), sqrt(2*(1+sin(t)))*sin(s)]
(which is 2D to 2D) maps the picture onto the desired Lambert projection.  This is shown in plot p2. (Well, actually in plot p2 I have embedded the result into 3D in order to later form a composite plot, but the embedded plot is inherently 2D.)

I think it is possible, although I haven't done this myself since it is late at night, to go from the 2D picture to the 3D sphere, rotate the sphere as needed, and then apply the Lambert transform, and figure out the overall map from 2D to 2D.

>  restart; >  p1 := plot3d([cos(t)*cos(s), cos(t)*sin(s), sin(t)],     s=0..2*Pi, t=-Pi/2..Pi/2,     image="earth.jpg",     labels=[x,y,z], size=[300,300]): >  p2 := plot3d([sqrt(2*(1+sin(t)))*cos(s), sqrt(2*(1+sin(t)))*sin(s), -2],     s=0..2*Pi, t=-Pi/2..Pi/2,     image="earth.jpg", size=[300,300]): >  plots:-display([p1,p2],     scaling=constrained, orientation=[145,72,0],     axes=none);

Download lambert-projection.mw

You say "after solving..." but you haven't asked Maple to solve the equation.  The solution u:=... that you have entered manually is incorrect.  If you insist on entering your own solution, then apply the quadratic formula with a little bit more care.  Alternatively, ask Maple to produce the solution for you.

See if this is useful.

 Download mw2.mw 

What you are asking is a job for your operating system, not Maple.

Install fswatch which is available for Linux, Mac, and Windows. Run it in a terminal and watch its notifications.  Then do what needs to be done in Maple.

I don't know how to do that with the tools provided in the VariationalCalculus package, but this can be done with bare hands with the help of the extended diff operator from the Physics package.  Have a look at this worksheet:   mw.mw

 

Download mw.mw

It this, and similar problems, look for the equilibria first.  That tells you where the interesting things happen.  In your case, the equilibria are at y=0 and y=1/9=0.111.  Therefore it makes sense to limit the y direction to the range 0 to 0.2.  Negative y values are not relevant to this problem since populations cannot be negative.  Negative time values are also of no interest since we are interested in predicting the future based on what we have now.  These considerations limit the plotting region to t > 0 and 0 < y < 0.2.  The upper limit for t is determined by trial-and-error through the expectation that the solution should converge to the stable equilibrium point.

I have modified your worksheet to account for those comments.  Here is what we get:

Worksheet: mw.mw

In Linux, one normally starts Maple through the command

maple -x &

In cases where the java memory is apt to be exhausted (often due to an extensive animation), one can request an increased  Java heap by changing that to

maple -x -j 4096 &

I think that option is undocumented.   The default value of the Java heap is 512.

That's how it works in Linux.  There may be similar hooks in other operating systems.

If you want insulated ends, you need to supply Neumann, not Dirichlet, conditions.  I have made that change in the worksheet below.

>  restart; >  pde := diff(u(x,t),t) = a^2*diff(u(x,t),x,x); >  ibc := D[1](u)(0,t)=0,  D[1](u)(L,t)=0, u(x,0)=x*(L-x); >  %pdsolve(pde, {ibc}, numeric): eval(%, {a=1, L=1}):  # set the parameter values pdsol := value(%); >  pdsol:-animate(t=0..0.3, frames=40);

Download mw.mw

Added later

Accuracy issues

Looking a bit more closely to the solution presented above, I see some disconcerting issues.  We know that insulated boundary conditions imply that the heat equation preserves the average of the solution.  Since the initial condition is x*(L-x), the average is L^2/6.  Therefore the solution should approach the constant value L^2/6.  In our case L=1, therefore that constant should be 1/6 which is approximately 0.1667.  In the animation shown above the solution converges to something like 0.12 which is very far from the target.

We may increase pdsolve's accuracy by reducing the spacestep parameter whose default value is L/20.  To get a limit which is reasonably close to 0.1667, we need to take spacestep=L/1000 (see the demo below) which is a very inefficient way of solving a PDE.  There may be options to pdsolve to obtain the solution more efficiently but I can't tell.  Perhaps someone else who is more knowledgeable about pdsolve's options can chime in.

Something else that I must point out is that the wavy oscillations in the solution curves produced by the default spacestep is a numerical artifact.  Note that those oscillations  are not present in the more accurate solutions.

Download mw2.mw

Download mw.mw

Note added later

The solutions shown above are correct but they can be simplified (by hand, not Maple!)  See the worksheet below for the simplification.

Download mw2.mw

A field plot for a nonautonomous system doesn't make much sense.  It is like asking for a picture of a stormy sea.  You may take a snapshot at one moment, but that picture will change at the next moment—there is no one picture.

It is possible, however, to convert a nonautonomous ODE to an autonomous one at the expense of introducing an extra dimension.  Specifically, one introduces a new unknown, let's say z(t), through the equation dz/dt=1, z(0)=0, which is equivalent to z(t)=t.  Then we replace t by z(t) on the right-hand side of the equation,.

Your equation is of the second order.  We convert it to a first order system, and then append the equation dz/dt=1 and thus obtain an autonomous system of 3 equations and then produce the corresponding 3D field plot.

See the details in the worksheet.

Download mw.mw

PS: Your title says "solve" but the code fragment that you have shown attempts to do a field plot which is something quite different.  In what I have written above, I have addressed the field plot issue.  If yiou are interested in a solution, then you will need to supply an initial condition.

I will show the solution of your homework problem's first part which can be done in Maple.  The rest, which calls for a straightforward mathematical argument, will have to be done by hand.  I have provided an outline below.

restart;
eq := diff(z(t), t) = (a + b*I)*z(t) + G(abs(z(t))^2)*z(t);
Eq := eval(eq, z(t)=r(t)*exp(theta(t)*I));
de1 := evalc(Re(Eq));
de2 := evalc(Im(Eq));

Then:

1. Solve {de1, de2} as an algebraic (not differential) system of two linear equations in the two unknowns {dr/dt, dtheta/dt}.  (This step can be done in Maple.)
2. Apply the assumptions on a, b, and G to conclude that the origin is the only the equilibrium.
3. Show that the origin repels orbits.
4. Show that the origin attracts far-away orbits.
5. Apply Poincare-Bendixon to conclude that there exists a limit cycle.