Here is an example, that does not require loading any packages.  There are other ways to do this, including pointplot, and scatterplot.   ydata is your "array" of points.  In your case it will have 25 elements. 

restart:

ydata:=[4.019,3.986,3.883,3.640,3.344,3.113,2.940,2.731,2.464,2.202,2.016,1.843];

xdata:=[seq(i,i=0..nops(ydata))];

pair:=(x,y)->[x,y]:

P:=zip(pair,xdata,ydata):

plot(P);

You can do it using the "display" command after a call to odeplot.   Here is an example taken from the help pages for the odeplot command.  Your post did not come through as intended. 

> restart:
> with(plots):
> p:= dsolve({D(y)(x) = y(x), y(0)=1}, type=numeric, range=-5..2):
> q:= dsolve({D(y)(x) = y(x)+3, y(0)=1}, type=numeric, range=-5..2):
> pp:=odeplot(p,color=blue):
> qq:=odeplot(q):
> display(pp,qq);
 

You can do it like this, but the answer returned is 0:

solve( {I1,I2,I3,I4}, M );

For more information regarding this, take a look at the help files for the solve command.

?solve

Look at the help files for "intersectplot" or type ?intersectplot and hit the enter key.

More details needed.  Are you looking for a 3d plot of x,y,z with a "box" around it like this?

with(plots):

implicitplot3d(x^3+y^3+z^3+1=(x+y+z+1)^3, x=-2..2, y=-2..2, z=-2..2, axes=boxed);

or a plot of a "box" in 3d like this?

with(geom3d):

cube(Ic,point(o,3,3,3),1): draw(Ic,axes=boxed);

Here is one way:

ans:=solve(x^4=16,x);

seq({m=ans[i],m=ans[i],m=ans[i],m=ans[i]},i=1..4);

{m = -2}, {m = 2}, {m = 2*I}, {m = -2*I}

I exported your worksheet as a text file and then looked at it.  There were issues with using code as text.  When inserting text, use the text button.  As it is, it loads without errors and runs without errors.

Download 221_16570_1.mw

Look at the help files for the plot command.  Type, or cut and paste the code below, and hit the enter key, or open the help browser and search for plot.

?plot

This is how you would define a function of x, using your first example.  Notice that I put a multiplication sign after the 9.

f:=x->x^3+3^2-9*x+5;

so after doing the above, if you enter;

f(3)

Maple will answer with

14

The help files have examples exactly like what you are trying to do.  Start there, and see how far you get.

restart:

v:=<2,4,5,8,9> ;

map(x->1/4*Pi*x^2, v);

evalf(%);

In Maple, you use exp(some stuff) instead of e^(some stuff), so the fit command would change to look like this:

fit = [a*exp(b*x), x]

?exp for more help on this

The Statistics package has commands that support plotting of data and curvefitting together.  See the help files for more info.

?ScatterPlot

restart:  with(plots):

with(Statistics):

 X := Vector([1, 2, 3, 4, 5, 6], datatype=float):
 
Y := Vector([2, 3, 4.8, 10.2, 15.6, 30.9], datatype=float):
 
R := ScatterPlot(X, Y, fit = [a*x^b, x], thickness = 3):
 
display(R);

The answer to "is the Fit command utilizing least squares? " is yes.

To display both, you can do it like this.  Just add the lines below to what you already have.

with(plots):

data:=plot(X,Y,style=point,symbol=circle,symbolsize=20,color=blue):

eqplot:=plot(eq,x=0..6):

display({data, eqplot});

Here is one way using the Statistics package.

restart:

with(Statistics):

X := Vector([1, 2, 3, 4, 5, 6], datatype=float):

Y := Vector([2, 3, 4.8, 10.2, 15.6, 30.9], datatype=float):

Fit(a*x^b, X, Y, x);

restart:

v:=t->diff(x(t),t);

PE:= 3*x(t)^2 + 4*x(t)^4;

KE:= 2*v(t)^2;

TE:=PE+KE;

plot(TE,x);

Does this work for you when you use your equations?

Review your original definition of J.  There is a section where you have a floating point value entered:  -3.0-3*n4*n10*n8*n2+3*n4*n10*n7*n1

If you get rid of the floating point value (-3.0) and replace it with a fixed point value (-3), then you get a solution fairly quickly.  Hope this helps.  Again, review your original input and make sure you entered it correctly.