select(t -> has(t,a) and has(t,b) and has(t,c), L);

Whether it's Fortran or C++ is not important, the question may be what method is being used.  Without seeing your system, it may be hard to say what's happening.  What does Maple's solution look like for, say, t =  2939.9?  Do you seem to be getting very large values or very large derivatives?  That could indicate that Maple is right that there is a singularity nearby.  You might try using a different method, or changing some of the options relating to error control (see ?dsolve,Error_Control).

RootOf(expr) (where expr is an expression containing the variable _Z) stands for a root of expr, i.e. any value of _Z that makes expr = 0.  In this case, those values are sqrt(58) and -sqrt(58).  In a situation that might involve different roots of the same equation, the label is useful in telling you that those with the same label must be the same. 

You might find the allvalues command useful here.

> allvalues(%);

[[x = ln(4500*sqrt(58)*(1/29)), y = 2*ln((1/3)*sqrt(58))]], [[x = ln(-4500*sqrt(58)*(1/29)), y = 2*ln(-(1/3)*sqrt(58))]]

> S := fsolve((x-h)^2 + (eq-k)^2 - r^2, x);
  if type([S],list(numeric)) and nops([S]) >= 1 then 
    this
  else
    that
  end if;

I assume what you have is a Matrix.

> with(LinearAlgebra):
  A:= RandomMatrix(8, 3);
  B:= SubMatrix(A, [seq(2*i-1,i=1..4)],1..3);

I assume l is the variable you want to solve for. 
In general, it is very unlikely that there is a closed-form solution. 

You could find numerical solutions (e.g. with fsolve) for numerical values of the parameters, or you might find series solutions in some situations.

If you don't want any solution curves, use dfieldplot instead.

> with(plots): with(DEtools):
 hyperb:= implicitplot(y^2 - x^2 - 1, x = -4 .. 4, y = -3 .. 3, colour=blue):
 phasep:= dfieldplot([(D(x))(t) = -.1+x(t)^2-x(t)*y(t), (D(y))(t) = y(t)^2-x(t)^2-1], [x(t), y(t)],  
     t=-1..1, x = -4 .. 4, y = -3 .. 3):
 display([hyperb,phasep]);

There are no such A and B if f and g are not coprime as polynomials in x.  If they are coprime, you can use Euclid's algorithm to find them.

The numeric PDE solver pdsolve(..., numeric) is rather limited.  As the error message says, it can only handle two independent variables.  Yours involves three: x, y and t.  Fortunately, the only place derivatives with respect to x or y appear is in the form diff(i(x, y, t), x)+diff(i(x, y, t), y).  Thus if you made a change of variables, say x = y + u, you would obtain a system with no derivatives with respect to u.  For each fixed u, you would then have a system with two independent variables
  

> EDP := [diff(s(y,t),t) = -.45*i(y,t)*s(y,t), diff(i(y,t),t) = .45*i(y,t)*s(y,t)-.2*i(y,t)-.7*diff(i(y,t),y)]
  IBC:= {s(-u, t) = 1, s(y, 0) = f(y+u, y), i( y, 0) = g(y+u)}; 

which pdsolve(..., numeric) should be able to handle.  For example:

> M:=pdsolve(EDP, eval(IBC, u=20), [s,i], numeric, time=t, range=-20..100);

You are defining phi[1](z) rather than defining the function phi[1].  What you probably meant is

> phi[1] := z -> A[1]*z^(lambda+I*epsilon)+B[1]*z^(lambda-I*epsilon);  

The way you have it, phi[1](z) is defined (technically, it's an entry in the remember table of phi[1]), but this would not define phi[1] of anything else.
The 2-D math parser is asking you whether you really meant what you wrote (the remember table entry) or the function definition as I wrote it.

In general, optimizing over an infinite set of integers is a difficult problem.  If you're willing to restrict to a finite set (small enough to compute all the values), it's easy.
 

> V:= [seq(f, k=1..10000)]:
  maxv:= max(V);  # this is the maximum value
  member(maxv,V,'kmax'):
  kmax;  # this is the position of the maximum value in the list

                         maxv := 0.9455376851

                                  100

 

Here's how I'd do it.   I'm using the parameter values suggested in the Wikipedia article on Parrondo's paradox.

Game A: payoff is +1 with probability 0.495, -1 with probability 0.505.

Game B: when current fortune is divisible by 3, payoff is +1 with probability .095, -1 with probability .905. 
Otherwise payoff is +1 with probability .745, -1 with probability .255.

Game C: play game A when turn number is divisible by 3, otherwise game B.

> with(Statistics):
  Coin1:= Sample(Bernoulli(.495), 20000):
  Coin2:= Sample(Bernoulli(.095), 20000):
  Coin3:= Sample(Bernoulli(.745), 20000):
  A[0]:= 0: B[0]:= 0: C[0]:= 0:
  for t from 1 to 20000 do
    d1:= 2*round(Coin1[t])-1;
    d2:= 2*round(Coin2[t])-1;
    d3:= 2*round(Coin3[t])-1;
    A[t]:= A[t-1] + d1;
    if B[t-1] mod 3 = 0 then
      B[t]:= B[t-1] + d2
    else
      B[t]:= B[t-1] + d3
    end if;
    if t mod 3 = 0 then
      C[t]:= C[t-1] + d1
    elif C[t-1] mod 3 = 0 then
      C[t]:= C[t-1] + d2
    else 
      C[t]:= C[t-1] + d3
    end if
  end do:
 plot([[seq([i,A[i]],i=0..20000)],[seq([i,B[i]],i=0..20000)],[seq([i,C[i]],i=0..20000)]],
   symbol=point,style=point,colour=[red,blue,green]);

Game A results are in red, game B in blue, game C in green.

That's an interesting expression you have there...  I hope I copied it right.

> Digits:= 40:
 Q:=.579920000000000000000000000000e-2*((-.9991*
(.38*sigma+10*l*coth(.862187887984549593047316871293*l))/sinh(l)
/(.7875e-2*cosh(.137812112015450406952683128707*l)+coth(.862187887984549593047316871293*l)
*sinh(.137812112015450406952683128707*l))*(coth(.862187887984549593047316871293*l)-coth(l))+9.9910
/sinh(.862187887984549593047316871293*l))*l*cosh(.862187887984549593047316871293*(l^2+sigma)^(1/2))
-(-.9991*(.38*sigma+10*l*coth(.862187887984549593047316871293*l))/sinh(l)
/(.7875e-2*cosh(.137812112015450406952683128707*l)+coth(.862187887984549593047316871293*l)
*sinh(.137812112015450406952683128707*l))*(coth(.862187887984549593047316871293*l)-coth(l))+9.9910
/sinh(.862187887984549593047316871293*l))*l*cosh(.862187887984549593047316871293*l)-sigma)
/sinh(.862187887984549593047316871293*(l^2+sigma)^(1/2))/sigma*cosh(.862187887984549593047316871293
*(l^2+sigma)^(1/2))*(l^2+sigma)^(1/2)+.579920000000000000000000000000e-2/sigma
*(-.9991*(.38*sigma+10*l*coth(.862187887984549593047316871293*l))
/sinh(l)/(.7875e-2*cosh(.137812112015450406952683128707*l)
+coth(.862187887984549593047316871293*l)*sinh(.137812112015450406952683128707*l))
*(coth(.862187887984549593047316871293*l)-coth(l))+9.9910
/sinh(.862187887984549593047316871293*l))*l*(-sinh(.862187887984549593047316871293*(l^2+sigma)^(1/2))
*(l^2+sigma)^(1/2)+sinh(.862187887984549593047316871293*l)*l)+.181225000000000000000000000001e-3
*(6*(-8.742125*l*(.38*sigma+10*l*coth(.862187887984549593047316871293*l))/sinh(l)
/(.7875e-2*cosh(.137812112015450406952683128707*l)+coth(.862187887984549593047316871293*l)
*sinh(.137812112015450406952683128707*l))*(-7/6+1/6*cosh(4*l)+cosh(l)^2
-1/3*sinh(l)*sinh(3*l))*cosh(.137812112015450406952683128707*l)+(sinh(l)-1/3*sinh(3*l))*sigma)
*(1.00000000000000000000000000000*l^2+.563063063063063063063063063063e-2*sigma)^(1/2)
*cosh(.137812112015450406952683128707*(1.00000000000000000000000000000*l^2
+.563063063063063063063063063063e-2*sigma)^(1/2))+52.452750*(-7/6+1/6*cosh(4*l)
+cosh(l)^2-1/3*sinh(l)*sinh(3*l))*(.38*sigma+10*l*coth(.862187887984549593047316871293*l))
/sinh(l)/(.7875e-2*cosh(.137812112015450406952683128707*l)+coth(.862187887984549593047316871293*l)
*sinh(.137812112015450406952683128707*l))*((1.00000000000000000000000000000*l^2
+.563063063063063063063063063063e-2*sigma)^(1/2)+l*sinh(.137812112015450406952683128707
*(1.00000000000000000000000000000*l^2+.563063063063063063063063063063e-2*sigma)^
(1/2))*sinh(.137812112015450406952683128707*l))*l)/(-1/2+cosh(l)^2-sinh(l)*cosh(l))/sigma/sinh(l)^3
/sinh(.137812112015450406952683128707*(1.00000000000000000000000000000*l^2+
.563063063063063063063063063063e-2*sigma)^(1/2))/(cosh(l)^2-1/2+sinh(l)*cosh(l))+.9991*
(.38*sigma+10*l*coth(.862187887984549593047316871293*l))*coth(.862187887984549593047316871293*l)
/(.7875e-2*coth(.137812112015450406952683128707*l)
+coth(.862187887984549593047316871293*l))-9.9910*l*coth(.862187887984549593047316871293*l);

> QS:= series(Q, sigma = 0, 2);

It's not easy to see, but this is of the form A*sigma^(-1) + B + O(sigma).
 

> A:= coeff(QS, sigma, -1);

Looks complicated, but...

> simplify(convert(A, exp)) assuming l > 0;
  

0

The same simplification works assuming l < 0 (but for l = 0 you run into a division by 0). 

So the sigma^(-1) term isn't really there, and the limit will be

> B:= coeff(QS, sigma, 0);

> V:= simplify(convert(B, exp)) assuming l > 0;

That very complicated expression in l took a long time to simplify, and was still very complicated 

> plot(V, l = 0 .. 5);

If the surface is defined by an implicit equation, just substitute the parametric form of the space curve into the equation for the surface and solve.

For example: where does the space curve (x = t, y = t^2-1, z = t^3-1) intersect the surface x*y*z = 1?

> curveq:= [x=t, y=t^2-1, z=t^3-1];
  surfeq:= x*y*z = 1;
  tvalues:= fsolve(subs(curveq, surfeq));

    tvalues := -1.157143622, 1.291980321

Since the equation in this case is polynomial, fsolve returns all the real solutions for t.  To visualize:

> with(plots):
   display([implicitplot3d(surfeq,x=-3..3,y=-3..3,z=-3..3,grid=[30,30,30],
      style=patchnogrid),
    spacecurve([t,t^2-1,t^3-1],t=-2..2,colour=black),
    pointplot3d(map(v -> eval([t,t^2-1,t^3-1],t=v),[tvalues]),symbol=box,
       colour=black)],
    view=[-3..3,-3..3,-3..3], axes=box, scaling=constrained, 
        orientation=[-100,70]);

As your examples show, it's generally a bad idea to use floats to index a table, because of this issue (and others involving roundoff error). 
But you might use this procedure to remove trailing zeros from a float:

> RemoveTrailingZeros:= proc(x) local S; S:= sprintf("%s0.%dg","%",Digits); parse(sprintf(S,x)); end proc;

> RemoveTrailingZeros(0.20);

         0.2