If you have a some guesses as to the general form of the curve and you are looking for a curve that approximates the data but does not necessarily pass through every point, use Statistics:-NonlinearFit or Statistics:-LinearFit. It's hard to say more without seeing your data.

This is a very difficult problem, and in Maple several intricate steps are required to solve it. I am amazed that it was assigned as a homework problem. Here's an overview of the steps:

1. Get a parameterized numeric dsolve solution with the initial velocity and and air density as parameters.
2. For any given density and initial velocity, extract the numeric x- and y-coordinate functions from the dsolve solution.
3. Use fsolve on the x-coordinate function to determine the time required for the ball to get to the home-run distance.
4. Put that time into the y-coordinate function to determine whether the ball clears the home-run fence.
5. Repeat steps 2-4 with different initial velocities to determine the minimal initial velocity needed to clear the home-run fence. (Unfortunately, fsolve is too fussy to perform this step. I had to resort to a crude bisection root-finding method.)
6. Repeat steps 2-5 with a different air density.

By Boyle's Law, pressure is proportional to density, so 10% lower air pressure corresponds to 10% lower air density.

Download homerun.mw

Try using a restart command. I am able to execute your Explore command without any problem.

Use an epsilon (I used 1 below) to change all your strict inequalities to nonstrict. Then pass them to Optimization:-LPSolve with any objective function (I used 0 below). 

e:= 1:
Optimization:-LPSolve(
     0,
     {x1 >= 0+e, x1+x2 >= 300+e, x2 <= 1000-e,
     x1+x2 <= 700-e, x2 >= 0+e, x1 <= 1000-e}
);

I think that my approach to making a color bar legend is sufficienty different from those found in the links that Alejandro provided that it is worth posting; however, it suffers the same drawback as all the others---that Maple's array plotting function allocates equal space to the color bar and to the main plot. The major difference is how I use plottools:-transform to extract the z-range from the plot.

First, make your plot and assign it to a variable. It can be any 3D plot produced with shading= zhue; it does not need to be produced with plot3d.

plot3d(x+y, x= 0..1, y= 0..1, shading= zhue);
P:= %:

(Plot omitted for brevity.) The following appliable module extracts the minimum and maximum z-coordinates from a 3D plot.

minMax:= module()
local min, Max;
export
     filter:= proc(x,y,z)
          if z::numeric then
               if z < min then min:= z end if;
               if z > Max then Max:= z end if
          end if;
          [x,y,z]
     end proc,
     ModuleApply:= proc(P::specfunc(anything, PLOT3D))
          min:= infinity;  Max:= -infinity;
          plottools:-transform(filter)(P);
          `..`(min,Max)
     end proc
;
end module;

Now use the plot's data to make the color bar:

plots:-tubeplot(
     [0,0,z], z= minMax(P),
     shading= zhue, scaling= constrained, radius= .2, axes= frame,
     tubepoints= 3, style= patchnogrid, orientation= [90,90],
     axis[1,2]= [tickmarks= 0], lightmodel= none
);
colorbar:= %:

Now use an array plot to put the two next to each other:

plots:-display(< P | colorbar >);

(Unfortunately, MaplePrimes cannot handle array plots.)

The following will produce a very smooth curve that fits the data rather closely.

m:= min(y_data): M:= max(y_data):
Y:= [seq(m..M, (M-m)/100)]:
X:= CurveFitting:-ArrayInterpolation(y_data, x_data, Y, degree= 3, method= spline):
plot(zip(`[]`, Y, X));

Since you said that you don't quite understand the math behind it, here it is in more detail.

F:= (x,y)-> x*(x+y)*exp(y-x):

We first note that the function is differentiable for all x and y.

Compute all first partial derivatives.
dFdx:= simplify(D[1](F)(x,y));

dFdy:= simplify(D[2](F)(x,y));

Set first partial derivatives to 0 and solve simultaneously. In the language of analytic geometry, this corresponds to finding where the tangent plane is horizontal. In the language of vector calculus, we would say that these points are where the gradient is 0. These points are called critical points.

CritPts:= [solve({dFdx=0, dFdy=0}, {x,y})];

Now we need to classify the critical points as minima, maxima, or saddle points. Algebraically, the easiest way to do this is to use the second derivatives. Since we have a computer algebra system, that's the way that we'll use.

Compute the matrix of second partial derivatives, which is called the Hessian.

Hess:= Matrix(2, (i,j)-> simplify(D[i,j](F)(x,y)));

We evaluate the Hessian at each critical point and then compute the eigenvalues. If all the eigenvalues are positive, the point is a minimum. If all the eigenvalues are all negative, it's a maximum. If some eigenvalues are positive, some are negative, and none are zero, then it's a saddle point. If any eigenvalues are zero, the test is inconclusive.

for cp in CritPts do convert(evalf(LinearAlgebra:-Eigenvalues(eval(Hess, cp))), list) od;
                  [2.414213562, -0.414213562]
                  [0.3543123712, 0.0516934784]

So (0,0) is a saddle point and (1/2, -3/2) is a minimum.

(Man, I'm sick and tired of posting whole worksheets by cutting and pasting each paragraph separately! Please fix this so that I can upload whole worksheets! Sometimes it works for me, sometimes not.)

The command Maplets:-Elements:-TextField has an option tooltip= "...". That's what you want. See ?Maplets,Elements,TextField

I wrote you a module (actually a module-valued proc) to help with the process of building an animation:

Frames:= proc()
    module()
    local frame:= 0, Frames:= table();
    export
         Add:= proc(P::seq(specfunc(anything, {PLOT,PLOT3D,ANIMATE})))
              frame:= frame+1;  
              Frames[frame]:= `if`(nargs=1, P, plots:-display(P))
         end proc,
         Animate:= ()-> plots:-display(convert(Frames, list), insequence),
         New:= proc()  frame:= 0;  Frames:= table()  end proc
     ;
     end module
end proc:

At the beginning of your worksheet, or somewhere before the first plot that you want included in the animation, include the line

GamesAnim:= Frames():

(The GamesAnim can be any name that you like.) Every time that you generate a frame (any plot that you want in the animation), use the module export Add:

GamesAnim:-Add(P1,P2);

When you want the animation, use

GamesAnim:-Animate();

If you want to start a new animation using the same kernel memory, use

GamesAnim:-New();

Here is your worksheet with these features added:

posterior_graphs_(animated)_1D.mw

Attempting to use solve to isolate the derivative provides a clue. For some inexplicable reason, the pieces of the piecewise are re-expressed as one-element lists rather than as simple expressions. After converting them back to expressions, the DEplot works fine with arrows.

Download DEplot-bug-solved.mw

Use the define command:

define(L, orderless, multilinear);
L(2*y1(x)+3*y2(x), 5*z1(x)+7*z2(x));



Many options exist for this command.

The reason is that everything is optimized for hardware floating-point double precision. That's what you get when you set Digits to 15 or fewer or set the Matrix's datatype to float[8]. This is not Maple's doing. Your processor is optimized for it, the compilers are optimized for it, etc.

If M is your Matrix, then

Statistics:-StandardDeviation(M[1, ..]);
Statistics:-Histogram(M[1, ..]);

Here's an example of the type of recursion that Acer referred to in his comment. Once again, I apply it to a Cartesian product.

CartProdRecursive:= proc(L::list(list))
     local a,b;
     if nops(L) < 2 then [L]
     else thisproc([[seq(seq([op(a),b], a= L[1]), b= L[2])], L[3..-1][]])[]
     end if
end proc:
     
CartProdRecursive([[a,b], [1,2], [c,d]]);

This is a very inefficient way to generate a Cartesian product. Perhaps Acer has in mind another example where this recursion technique is competitive with other techniques.

Instead of "Modulus of Elasticity for Element i", use cat("Modulus of Elasticity for Element ", i).