Replace the last line by the line

solve({seq(lhs(U)[i] = rhs(U)[i], i = 1 .. 2)}, {Gx, Gy})

Solve command does not solve vector equations. Index  i  not equals to 3, because the last coordinate of lhs(U) does not contain  Gx  and  Gy .

Of course, we assume that  lhs(U)[3]=0 , otherwise there are no solutions.

Assignment   i:='i'  is incorrect, since then  Term1:=i-1  has not a specific value, and therefore the command  add  can not be executed.

algsubs command  is more powerful than applyrule command.

Example:

a:=x[2]+y-b+x[1]-c+x[5]+d+x[4]+x[3]+expand(sum(x[i], i=1..5)^2);
applyrule(sum(x[i], i=1..5)=Sum(x[i], i=1..5),a);
algsubs(sum(x[i], i=1..5)=Sum(x[i], i=1..5), a);

u:=sin(w*t-theta)+sin(w*t-theta-2*Pi/3)-sin(w*t-theta+2*Pi/3);
algsubs(w*t-theta=x, u);
expand(%);
simplify(applyrule(a::algebraic*sin(x)+b::algebraic*cos(x)=sqrt(a^2+b^2)*sin(x+arctan(b,a)),%));
subs(x=w*t-theta, %);

is here  http://www.mapleprimes.com/questions/145061-Based-On-2-Sets-Of-Vectors-In-R4-How 

Read help  ?LinearAlgebra[LinearSolve]

One vector equation  x1*v1+x2*v2+x3*v3=x4*w1+x5*w2+x6*w3  is equivalent to the homogeneous system of four scalar equations with six unknowns from  x1  to  x6 . Since the rank of the leading matrix is greater  or equal to 3, the dimension of the solution space will be equal to 2 or 3. Linear systems in Maple are solved by  LinearAlgebra[LinearSolve]  command.

without changing the code of the procedure:

a:=1: b:=2:

w(a,b); 

         99

Your problem can be easily solved without Maple by using an iteration method for solving equation  x=f(x)  (problem on fixed point of mapping  f ). Here  f:=x->(x+2)/(x+3) .  See  http://en.wikipedia.org/wiki/Fixed_point_(mathematics) 

All conditions for the convergence of your iterations  a_n = (a_(n-1) +2)/(a_(n-1) + 3),  a_1 =2   performed:

1) Function  f  maps a closed segment  [0, 2]  into itself.

2) On this segment  [0, 2]  mapping  f  is a contraction, since  |f'(x)|=1/(x+3)^2<=1/9

So your iteration converges to the unique root of the equation   x=(x+2)/(x+3)  on the segment  [0, 2], ie to the number  sqrt(3)-1 .

By default, the  Vector  command means a column vector. Even shorter, you can specify it if you use angle brackets instead Vector  command:

Nr:=5:

p:=2*<(5-i $ i=0..Nr)> + <(1 $ Nr+1)>;

Painting.mw

restart;

with(DEtools): with(plots): Nsols := 5; Ntstep := 10;

 k := 0; A := 0.37e-1; B := 0.2e-6;

ode0 := diff(U(t), t) = -(A+B*U(t))*U(t);

 ic[1] := U(365*k) = 1000;

sol[0] := dsolve({ode0, ic[1]}, U(t), range = 365*k .. 365*k+365, numeric);

p[0] := odeplot(sol[0], [[t, U(t)]], t = 365*k .. 365*k+365, colour = blue, legend='p[0]'):

display(p[0]);

ode1 := diff(U(t), t) = -(A+r(t)+B*U(t))*U(t);

R := RandomTools:-Generate(distribution(Normal(-0.2e-1, 0.4e-1)), makeproc = true);

r := proc (t) if not type(t, numeric) then return 'procname(args)' end if; R() end proc;

for i to Nsols do

sol[i] := dsolve({ode1, ic[1]}, numeric, known = r(t), method = classical[foreuler], stepsize = 0.1e-1):

p[i] := odeplot(sol[i], [[t, U(t)]], t = 365*k .. 365*k+365, colour = red)

 end do:

display(seq(p[i], i = 1 .. Nsols));

for j to Ntstep do

 for i to Nsols do

x[i] := proc (t) options operator, arrow; rhs(sol[i](t)[2]) end proc;

y[i] := x[i](j):

end do:

ymaxval[j] := max([seq(y[i], i = 1 .. Nsols)]):

end do:

S1:=[seq(ymaxval[j], j=1..Ntstep)];

for j to Ntstep do

 for i to Nsols do

x[i] := proc (t) options operator, arrow; rhs(sol[i](t)[2]) end proc;

y[i] := x[i](j):

end do:

yminval[j] := min([seq(y[i], i = 1 .. Nsols)]):

end do:

S2:=[seq(yminval[j], j=1..Ntstep)];

plotofmax:= plot([seq(j, j=1..Ntstep)], S1, color=red, legend="plotofmax"):

plotofmin:= plot([seq(j, j=1..Ntstep)], S2, color=black, legend="plotofmin"):

display(plotofmax, plotofmin, p[0], view=[1..Ntstep, min(S2, rhs(sol[0](10)[2]))..max(S1, rhs(sol[0](0)[2]))]);

Now all works:

restart;

with(DEtools): with(plots):

k := 0; A := 0.37e-1; B := 0.2e-6;

ode1 := diff(U(t), t) = -(A+r(t)+B*U(t))*U(t);

ic[1] := U(365*k) = 1000;

R := RandomTools:-Generate(distribution(Normal(-0.2e-1, 0.4e-1)), makeproc = true);

r := proc (t) if not type(t, numeric) then return 'procname(args)' end if; R() end proc;

for i to 5 do

sol[i] := dsolve({ode1, ic[1]}, numeric, known = r(t), method = classical[foreuler], stepsize = 0.1e-1):

p[i] := odeplot(sol[i], [[t, U(t)]], t = 365*k .. 365*k+365, colour = red) end do:

display(seq(p[i], i = 1 .. 5));

for i to 5 do x[i] := proc(t) options operator, arrow; rhs(sol[i](t)[2]) end proc:

y[i] := x[i](2) end do:

seq(y[i], i=1..5);

max(seq(y[i], i=1..5));

If you want a vector to be argument in your function, you can use the  unapply  command:

F:=unapply(6*x[1]^2 + 3*x[2]^2 - 1.5*x[3]^3, x::Vector):

Examples:

F(<1,2,3>);

H:=unapply(Student[VectorCalculus][Hessian](F(<X[1], X[2], X[3]>), [X[1], X[2], X[3]]), X::Vector):

H(<x,y,z>);

H(<1,2,3>);

To draw a phase curve, you must:

1) Set the value of the parameter  k .
2) Specify the initial condition.

If you interested of changing of the curve in response to changes  k , you can make animation with the parameter  k .

Example (for greater clarity, the ranges of  k  and  t  are reduced):

with(DEtools):
for k from 1 to 2 by 0.01 do
DE := diff(x(t), t) = x(t):
DF := diff(y(t), t) = k*y(t):
A[k]:=phaseportrait([DE, DF], [y, x], t = 0 .. 0.5, [[y(0)=1, x(0)=1]], color=black, linecolor = red):
od:
plots[display](seq(A[k], k=1..2, 0.01), insequence=true);

To construct the phase curves you need:

1) Increase the range for the variable t .

2) Set the initial conditions.

Example (plotting of two phase curves):

with(DEtools):
sys := diff(x(t), t) = y(t), diff(y(t), t) = x(t)*(1-x(t)*x(t))+y(t);
DEplot([sys], [x(t), y(t)], t = 0 .. 3, x = -3 .. 3, y = -3 .. 3, [[x(0)=1, y(0)=2], [x(0)=-2, y(0)=1]],color = black);