Be aware that the ?GraphTheory package does not permit graphs with self-loops.  A self-loop is implied if there are any non-zero elements on the diagonal of the adjacency graph.  For example:

GraphTheory:-Graph(Matrix(2, [[1,1],[1,0]]));
Error, (in GraphTheory:-Graph) an undirected graph cannot contain loops (loop
detected at vertex 1)

What version of Maple are you using, and is it on a 32 or 64 bit machine?  This works for me in in cmaple, the Standard GUI with both 1D and 2D input.

Your description is not precise.  The term 'function' in Maple specifically refers to an unevaluated function call, say f(x[1],3).   Is that what you mean?  If so, then where are the expressions of x[1] for which you want coefficients?  Inside the a, b, and c?   While a general description of the problem doesn't hurt, a specific example is usually easier to understand and allows one to test.
 

Maple tables are hash tables. 

A straightforward way is to use a table to store the data, then convert the table to the form you need.  For example:

temp := proc(xf, dx, dt)
local cnt, Data, x, t;
    Data := table();
    t := 0;
    x := 0;
    cnt := 1;
    Data[cnt] := [t,x];
    while x < xf do
        x := x+dx;
        t := t+dt;
        cnt := cnt+1;
        Data[cnt] := [t,x];
    end do;
    return [seq(Data[cnt], cnt=1..cnt)];  # the use of cnt as the index is a trifle tricky
end proc:

Data := temp(3, 0.1, 0.2);

plot(Data);

You can change the error tolerance using the abserr and relerr parameters.  You could also use a different integration method.

This is the chain-rule, using Leibnitz notation.  Because there is some ambiguity in the expression, Maple uses a different notation:

(**) z := x(y(t));
                                                              z := x(y(t))

(**) diff(z,t);
                                                                     /d      \
                                                          D(x)(y(t)) |-- y(t)|
                                                                     \dt     /

# Here we convert D(x)(y(t)) to the diff notation,
(**) map(convert,%,diff);
                                                    / d       \|          /d      \
                                                    |--- x(t1)||          |-- y(t)|
                                                    \dt1      /|t1 = y(t) \dt     /

You can use  multiline comments:

(*
   These lines won't be executed
*)

Alas, they are not available in 2D mode.

Why not just compute the cost for one year of usage?

Fluorescent = ($3.67)*4 + (13*W)*(4*hr/day)*(365*day)*(0.151e-3*dollar/W/hr) = $17.55
Incandescent = (3.29)/2  + (40*W)*(4*hr/day)*(365*day)*(0.151e-3*dollar/W/hr) = $10.46

It works fine here, in command line and in the Standard worksheet, both with 1D and 2D input.  Note that delta=1 does not result in alpha*t but rather alpha + t.

Consider trying fsolve.  For example,

f := exp(-x^2*sin(x)):
F := Int(f, x=0..b):
fsol := fsolve(F = 1/2, {b});
                                                       fsol := {b = 0.5166551554}
# check the result
evalf(eval(F, fsol));
                                                              0.5000000000

Here's a different approach.  It uses A[[i,j,k]] to select indices i,j,k from list (or set) A.

A:=[5,0,1,0,0]:
C:=[0.735,0,0.265,0,0]:
nonZeroIndices := [seq(`if`(A[i]=0,NULL,i), i=1..nops(A))];
                           nonZeroIndices := [1, 3]

A[nonZeroIndices];
                                    [5, 1]

C[nonZeroIndices];
                                [0.735, 0.265]

It would help if you didn't run separate equations together.  Punctuation can improve clarity.  I don't see how you came up y1=y1+/-exp(t*siga).

I found a site, www.mathcurve.com/fractals/gosper/gosper.shtml, that shows how to construct the base curve.  Here is an update that gives the precise curve. This can be cleaned up, but I took the easy route from what I'd already done.

recursive := proc(s::list, f::list, c::list, depth::posint)
local z,seg,segs;
    segs := seq(map(Transform, z, s, f), z in c);
    if depth = 1 then
        segs[];
    else
        seq(procname(seg[],c,depth-1), seg in segs)
    end if;
end proc:

Transform := proc(z,s,f)
local x,y,x2,y2,mag,theta;
    (x,y) := op(f-s);
    mag := sqrt(x^2+y^2);
    theta := arctan(y,x);
    x2 := z[1]*cos(theta) - z[2]*sin(theta);
    y2 := z[1]*sin(theta) + z[2]*cos(theta);
    [s[1] + mag*x2, s[2]+mag*y2];
end proc:

# Generate points of base curve
V := [seq([cos,sin](k*Pi/3), k=0..6)]:
E := [seq(V[k+1]-V[k], k=1..6)]:

P[1] := [0,0]:
P[2] := P[1]+E[5]+E[6]:
P[3] := P[2]+E[1]+E[6]:
P[4] := P[3]+E[3]+E[2]:
P[5] := P[4]+E[1]+E[2]:
P[6] := P[5]+E[5]+E[6]:
P[7] := P[6]+E[5]+E[6]:
P[8] := P[7]+E[5]+E[4]:

# Compute point (x,y) which is used to transform
# P[8] to (1,0); that is needed by Transform.

eqs := normal(Transform(P[8], P[1], [x,y]));
sol := solve({eqs[1]=1,eqs[2]=0}, {x,y});
xy := eval([x,y],sol);

# Apply Tranform to vertices
Q := map(Transform, [seq(P[k], k=1..8)], P[1], xy);

# Generate the Peano-Gosper curve
flowsnake := recursive([0,0],[1,0]
                       , [NULL
                          , [Q[1],Q[2]]
                          , [Q[3],Q[2]]
                          , [Q[4],Q[3]]
                          , [Q[4],Q[5]]
                          , [Q[5],Q[6]]
                          , [Q[6],Q[7]]
                          , [Q[8],Q[7]]
                         ]
                       , 3
                      ):

plot([flowsnake],axes=none,scaling=constrained);

Here's one approach.  I only approximated the shape, so you'll have to clean up

recursive := proc(s::list, f::list, c::list, depth::posint)
local z,seg,segs;
    segs := seq(map(Transform, z, s, f), z in c);
    if depth = 1 then
        segs[];
    else
        seq(procname(seg[],c,depth-1), seg in segs)
    end if;
end proc:

Transform := proc(z,s,f)
local x,y,x2,y2,mag,theta;
    (x,y) := op(f-s);
    mag := sqrt(x^2+y^2);
    theta := arctan(y,x);
    x2 := z[1]*cos(theta) - z[2]*sin(theta);
    y2 := z[1]*sin(theta) + z[2]*cos(theta);
    [s[1] + mag*x2, s[2]+mag*y2];
end proc:

# These are just a crude approximation

P0 := [0,0]:
P1 := [.3, -.125]:
P2 := [.60,  .125]:
P3 := [.2,   .250]:
P4 := [.15,  .60]:
P5 := [.575,  .5]:
P6 := [0.9,  .40]:
P7 := [1,     0]:

flowsnake := recursive([0,0],[1,0]
                       , [NULL
                          , [P0,P1]
                          , [P2,P1]
                          , [P3,P2]
                          , [P3,P4]
                          , [P4,P5]
                          , [P5,P6]
                          , [P7,P6]
                         ]
                       , 2 # increase to refine
                      );

plot([flowsnake],axes=none,scaling=constrained);