You usually need to solve for the same number of variables as there are equations. So, you need to pick the right 10 variables to solve for. It's clear from your "by hand" solution that you want to eliminate vx, vy, and the eight tgt variables. That's 10. So, if you do the following, you'll get exactly the same solutions as you got "by hand".

Sol:= solve(eqs, {vx, vy, seq([tgtX[k], tgtY[k]][], k= 1..4)}):
{tgtX[4]=eval(tgtX[4], Sol), tgtY[4]=eval(tgtY[4], Sol)};

You can give a list of points as the first argument of the regular plot command and by default you'll get an interpolated curve. It's just "connect the dots." So, in the simplest case, the command would be

plot(list of points);

No other arguments are needed.

Answering your Question directly, yes, there's an option that you can give to plots:-pointplot to achieve the same thing. That option is style= line.

What's needed is a better way to form Cartesian products. Using a nested seq or for loop can't be generalized to an arbitrary number of factors. Using combinat:-cartprod is awkward and too slow. The following is the best procedure from an extensive post on Cartesian product procedures by Joe Riel. The procedure is by Joe Riel, with some small modifications by me, mostly so that it would fit on one line.

CartProdSeq:= proc(L::seq(list), $)
local S,j;
    eval(subs(S= seq, foldl(S, [_i||(1..nargs)], seq(_i||j= L[j], j= nargs..1, -1))))
end proc:

Using it for your situation,

ex:= V[2] + a[5]*V[5] + a[6]*V[6] + a[7]*V[7] + V[9]:
v:= [a[5], a[6], a[7]]:
seq(eval(ex, v =~ k), k= [CartProdSeq(seq([0,x], x= v))]);

My favorite way is to use the trace command. For any (non-builtin) procedure P , you can give the command

trace(P):

before executing the procedure. Then when you execute code containing a call to P, you'll see when P is called, what its arguments are, what its return value is, and some of what happens inside P.

You can have trace turned on for any number of procedures at the same time.

Here are some somewhat elegant procedures to determine whether intervals intersect and what the intersection is. They work for any number of any type of real interval: open, closed, infinite, half-open, etc.

Download Intersect_Intervals.mw

With a few minutes (less than five) of pattern matching with the "naked eye", I came up with the following formula:

int(diff(P(n,x), x)*diff(P(m,x), x), x= -1..1) = `if`((m+n)::odd, 0, min(m,n)*(min(m,n)+1))

I verified the formula up to n=30, m=30. I suspect that the formula could be easily proved.

Here's how to create patterns for visual inpection:

interface(rtablesize= infinity):
M1:= Matrix(30, (m,n)-> int(diff(P(m,x), x)*diff(P(n,x), x), x= -1..1));

ColumnRelations:= proc(M::Matrix, {free::symbol:= :-_v})
uses LA= LinearAlgebra;
local k;
     evalc(map2(`.`, <seq(free[k], k= 1..op([1,2], M))>, LA:-NullSpace(M))) =~ 0
end proc:

ColumnRelations(M, free= v);

If the output in your worksheet can be regenerated quickly, then you should remove the output before saving it. You can use the menu selection Edit -> Remove output -> From worksheet.

Yes, you should consider dividing your worksheet. When the GUI gets bogged down, as often happens with lengthy output, there's no equivalent of the Stop button that will stop it. I've often had to kill the GUI, which kills the whole Maple process.

The worksheet below has the complete answers to your three questions, using a slightly different technique than Kitonum (very slightly different). I only post it because I was composing it at the same time that he was composing his.

Download zero-one_system.mw

Assuming that the input is floating point, the command that you want is frem(..., evalf(2*Pi)). If your input contains symbolic Pi that you'd like to keep in the output, let me know.

Here's a procedure that does what you want:

YourProc:= proc(
     fx::algebraic,
     `x= a..b`::And(name=range(realcons), satisfies(r-> op([2,1], r) < op([2,2], r))),
     n::posint,
     {functionoptions::list({symbol, name= anything}):= []},
     {pointoptions::list({symbol, name= anything}):= []}
)
local
     x:= op(1, `x= a..b`),
     a:= op([2,1], `x= a..b`),
     r:= op([2,2], `x= a..b`) - a,   #b-a
     Df:= unapply(diff(fx, x), x),
     k,   #index
     xs:= [seq(a+k*r/n, k= 0..n)]
;
     plots:-display(
          plot(fx, `x= a..b`, functionoptions[]),
          plot([seq([k, Df(k)], k= xs)], style= point, pointoptions[]),
          _rest
     )
end proc:

And an example of its use:

YourProc(
     x^2, x= -2..2, 5,
     functionoptions= [color= green],
     pointoptions= [symbol= cross, symbolsize= 24],
     axes= boxed
);

The last three arguments are optional. You could just as well use it as

YourProc(x^2, x= -2..2, 5);

Note that intervals in Maple are specified as 2..5 rather than [2,5]. However, if you insist on using [a,b], I can easily change the procedure to accomodate that.

As a workaround, you can use Maple Input (aka 1D input).

The direct cause of your problem is that you're not allowed to assign to a parameter, in this case M. You can make a local copy of M like this:

M1:= copy(M);

Then make all assignments to M1 and return M1.

There are also a few places where you used = that should've been :=.

You should begin your procedure with 

uses LinearAlgebra;

or

uses LA= LinearAlgebra;

Don't rely on the with statement being active in your procedure---sometimes it may be, sometimes it may not be.

Suppose that your dsolve result is assigned to ans:

ans:= dsolve(...);

Then do

u:= unapply(eval(u(t), ans), t);

If g=0, then f can be arbitrary and the jacobian will still be zero. That's why the fourth solution shows f unchanged.