Of course the second calculation  evalf(gg(2))=-0.3333333333  is incorrect. Probably Maple uses formal geometric series summation  b+b*q+b*q^2+..+b*q^n+..=b/(1-q)  which is correct only if  |q|<1 .

Indeed for the series  gg(2)=Sum(2^jj*2^jj, jj = 0 .. infinity)   we have  b=1 ,  q=4    and   1/(1-4)=-0.3333333333

You can search for a numerical solution using  numeric  option. To do this, all the parameters and initial conditions must be specified. With the resulting procedure, you can find solutions in some points, build plots, etc.

sys := {diff(x(t), t) = -k*x(t)*sqrt(x(t)^2+y(t)^2)/m, diff(y(t), t) = -k*y(t)*sqrt(x(t)^2+y(t)^2)/m - g}:

sol := dsolve({eval(op(sys), {g = 9.8, k = 1, m = 2}), x(0) = 0, y(0) = 0}, numeric);

plots[odeplot](sol, [[t, x(t)], [t, y(t)]], t = 0 .. 2, color = [red, blue], thickness = 3);

If you mean the 50 roots of this equation, then here are the first 50 positive roots:

seq(fsolve(0.5*lambda*tan(lambda)-1, lambda=`if`(n=0,0..Pi/2,-Pi/2+n*Pi..Pi/2+n*Pi)), n=0..49);

 Addition. You can easily plot these points with the following code (for clarity, plotted the first 10 roots in the form of small red circles):

L:=[seq(fsolve(0.5*lambda*tan(lambda)-1, lambda=`if`(n=0,0..Pi/2,-Pi/2+n*Pi..Pi/2+n*Pi)), n=0..49)]:

plots[display](plot(0.5*lambda*tan(lambda)-1, lambda=0..Pi/2+9*Pi+0.01, -10..10, color=blue, discont), plot(map(t->[t,0], L[1..10]), style=point, symbol=solidcircle, symbolsize=12, color=red, scaling=constrained));

Complex function of a complex argument may be regarded as a mapping from R^2 into R^2. On the left - the space of arguments,  on the right - the space of their images:

R:=evalc(1/(1-(x+I*y)));

R1:=op(1,R),coeff(R,I);

f:=unapply([R1],x,y);

F:=plottools[transform](f):

A:=plot([seq([convert(<cos(Pi/4*i),-sin(Pi/4*i);sin(Pi/4*i),cos(Pi/4*i)>.<t,0>,list)[], t=-0.8..0.8], i=0..7), seq([r*cos(t),r*sin(t),t=0..2*Pi], r=0.2..0.8,0.2)], color=[red,gold,brown,yellow,blue,green,violet,cyan], thickness=3):

plots[display](<A | F(A)>, scaling=constrained);

 Edited.

If we go to polar coordinates, Maple can easily find the solution in explicit form, even simpler than Mathematica does:

restart;

T:={x(t)=r(t)*cos(phi(t)), y(t)=r(t)*sin(phi(t))}:  # The transformation to polar coordinates

sys:={diff(x(t),t) = -k/m*x(t)*sqrt(x(t)^2+y(t)^2), diff(y(t),t) = -k/m*y(t)*sqrt(x(t)^2+y(t)^2)}:  # The original system

subs(csgn(r(t))=1, simplify(eval(sys,T)));  # The transformation of original system to polar coordinates

sys1:=solve(%,{diff(r(t),t),diff(phi(t),t)});  # The original system in polar coordinates in standard form

dsolve(sys1, {r(t),phi(t)});  # The general solution in polar coordinates

dsolve(sys1 union {r(0)=1,phi(0)=Pi/3},{r(t),phi(t)});  # The solutin with some initial conditions

eval(T, %);  # The same solution in Cartesian coordinates

In fact, the point moves along a straight line toward the origin.

ConvIntoIneq  procedure solves the problem:

ConvIntoIneq:=proc(Rel)

local t, R, x, y;

t:=indets(Rel)[];

R:=solve(Rel);

x, y:=op(1,R), op(2,R);

if x=-infinity then if type(y,realcons) then return cat(t, "<=", y) else

cat(t, "<", op(1,y)) fi else

if type(x,realcons) then return cat(t, ">=", x) else

cat(t, ">", op(1,x)) fi; fi;

end proc:

Examples of use:

ConvIntoIneq(not x<100), ConvIntoIneq(not x>100), ConvIntoIneq(not x>=100), ConvIntoIneq(not x<=100), ConvIntoIneq(x>100), ConvIntoIneq(x<100), ConvIntoIneq(x>=100) ;

 Addition. Unfortunately, if we want to get a "nice" output, the name of variable can appear on the right side of the inequality:

map(parse, [ConvIntoIneq(not x<100), ConvIntoIneq(not x>100), ConvIntoIneq(not x>=100), ConvIntoIneq(not x<=100), ConvIntoIneq(x>100), ConvIntoIneq(x<100), ConvIntoIneq(x>=100)])[] ;

I added 2 options into your code: axes=normal and orientation=[15,50]  and used user lighting :

 Addition. The plot was saved in Paint as png - file.

Another variant. Additionally I added  color=khaki ,  style=surface , numpoints=20000  and changed  user lighting  :

plot([[0,0.36],[5,0.38],[10,0.3],[15,0.18],[20,0.96],[25,0.18],[30,0.16],[35,0.96]], style=line, thickness=5, color=red, view=[0..35,0..1], axis[1] = [gridlines = [7, thickness = 1, subticks = false, color = grey]], axis[2] = [gridlines = [6, thickness = 1, subticks = false, color = grey]], tickmarks=[[0=1,5=5,10=10,15=50,20=250,25=1000,30=10000,35=50000], spacing(0.2,0)], size=[650,350], font=[COURIER,ROMAN,14]);

 Edited.

P := (L, a, b) ->plot(L, x = a .. b, color = [seq(`if`(i::odd, red, black), i = 1 .. nops(L))]):   # L is a list of functions

Example of use:

P([sin(x), cos(x), 2, 3-x], 0, 5);

First, you can find the indices of columns that must be deleted, and then delete them.

Example:

A:=<0,1,0,3; 0,-2,0,4>;

ListTools[SearchAll](0, convert(A[1], list));

LinearAlgebra[DeleteColumn](A, [%]);

Your syntax is wrong there. Below is the text of the procedure called  P  that solves the problem for any  Time  and  solution:

P:=proc(Time, solution)

local n;

n:=nops(Time);

piecewise( seq( op([t>=Time[i] and t<Time[i+1], solution[i]]), i = 1..n-1));

end proc:

Solution of above example:

P([0,1,2,3], [t^2, 1, 3-t]);

plot(%, t=0..3, scaling=constrained);

Edited.

Example:

F:=piecewise(t>=0 and t<1, t^2, t>=1 and t<2, 1, t>=2 and t<3, 3-t);

plot(F, t=0..3, scaling=constrained);

Addition:  outside the interval t=0..3  the function  F by default considered to be equal to 0 . 

Examples of plotting  4 basic trigonometric functions:

F1 := sin(x): F2 := cos(x): F3 := tan(x): F4 := cot(x): R := x = -Pi .. 2*Pi: C := [red, blue, green, brown]: c := scaling = constrained:

v1 := plot(F1, R, color = C[1]):

v2 := plot(F2, R, color = C[2]):

v3 := plot(F3, R, -3 .. 3, color = C[3], discont):

v4 := plot(F4, R, -3 .. 3, color = C[4], discont):

plot([F1, F2, F3, F4], R, -3 .. 3, color = C, discont, c);  # All together in one plot

 plots[display](< v1, v2; v3, v4 >, c, size = [300, 300]);  # Matrix of the plots

 plots[display](< v1 | v2 | v3 | v4 >, c);  # Vector-row of the plots

 plots[display](< v1, v2, v3, v4 >, c, size = [100, 100]);  # Vector-column of the plots

If you want the tickmarks are not lost while reducing the size of the plot, reduce the font size:

combine(y(x) = (1/3)*exp(x)+_C1/(exp(x))^2);

                                     y(x) = 1/3*exp(x)+_C1*exp(-2*x)

Example with other powers not exp:

combine(y(x) = 1/3/2^x+_C1/(3^x)^2);

                                      y(x) = 1/3*2^(-x)+_C1*3^(-2*x)