I'll answer this under the assumption that you'll want to call the optimized code with evalhf or compile it because that seems like the most likely scenario and there aren't many options for how things can be done. For one thing, procedures can't return Vectors, and they can't work with lists at all. To use an Array (Vector, Matrix) in evalhf, it must be passed in.

So, instead of making assignments to the entries of V, leave V undefined, and make equations that makeproc will turn into assignments.

restart:

L:= [V[1]= x+y,
         V[2]= x^2+y^2,
         V[3]= x^3+y^3,
         V[4]= x^2+y^2+x+y
       ]:

T:= codegen:-makeproc([codegen:-optimize(L, 'tryhard')], [x,y,V::'Vector']);

     T := proc(x, y, V::Vector)
     local t25, t26, t29, t30;
          t30 := x + y;
          t25 := y^2;
          t26 := x^2;
          t29 := t25 + t26;
          V[1] := t30;
          V[2] := t29;
          V[3] := t25*y + t26*x;
          V[4] := t29 + t30;
     end proc;

Note that V is just an undefined symbol throughout this process; I never assigned a Vector or anything else to it.

You essentially already gave the answer to your own question: Use either package Threads or package Grid. Without you providing more details, this is all the answer that you're likely to get:

If your procedures meet the stringent requirements of the Threads package (see ?index,threadsafe), use

Threads:-Add[tasksize= 1](f[k](x,y), k= 1..2);

Otherwise, use

`+`(Grid:-Seq(f[k], k= 1..2));

You said that member works for tables and lists, which is correct. However, ListTools:-SearchAll doesn't work for tables. So, in case you wanted to do this with a table, I wrote a procedure for it:

SearchAllTable:= proc(
     x::depends(`if`(sequence, list, anything)),
     T::table,
     {
          sequence::truefalse:= false,
          nolist::truefalse:= false,
          indexorder::{truefalse, procedure}:= false
     },
     $
)
local
     INDEXORDER:= `if`(
          indexorder::procedure,
          ':-indexorder'= indexorder,
          `if`(indexorder, ':-indexorder', [][])
     )
;
     indices(T, `if`(nolist, ':-nolist', [][]), INDEXORDER)
          [[ListTools:-SearchAll(`if`(sequence, x, [x]), [entries(T, INDEXORDER)])]]
end proc:

This has three optional arguments: nolist, indexorder, and sequence. The first two of these mean exactly the same thing as they mean as options to indices or entries (see ?indices). Since the entries of a table can be sequences (unlike the entries of a list), the option sequence means that the first argument is a list whose corresponding sequence is the search target.

Examples:

T:= table([k= seq(irem(k,4), k= 1..50)]):
SearchAllTable(3, T, nolist);

     3, 7, 11, 15, 19, 23, 27, 31, 39, 35, 47, 43

T[100]:= 4,5:
SearchAllTable([4,5], T, sequence);

     [100]

This is for the purpose of visual display only; it has no computational effect. Suppose that A is your Matrix. Then do

printf("%{c }a\n", A);

A preliminary experiment that I did suggests that you may be able to get a factor-of-3 time improvement by using codegen:-optimize(..., tryhard). This is an extreme form of computational simplification (which tends to make an expression visually atrocious and unintelligible---but that doesn't matter). I applied it to your expression g5 like this:

G5:= unapply(subs(Pi= evalf(Pi), g5), (t,x)):
G5O:= codegen:-optimize(G5, tryhard):
G5O:= subs(pow= `^`, eval(G5O)):
G5O:= subs(2^(1/2)= sqrt(2.), eval(G5O)):

The second of these steps takes about 30 seconds on my computer. The others are instantaneous.

Then I compared the times to evaluate g5 and G5O over a random sample of 512 pairs of floats uniformly distributed within your given ranges. I did this in separate sessions (i.e., with an intervening restart) to avoid remember table effects. The time for g5 was 12.52 s; the time for G5O was 4.34 s.

I've made no attempt to address the potential logic errors in your code that have been pointed out by Axel and Professor Hofri.

Setting parameter a=1 and making some assumptions, I get a result that includes a hypergeometric function.

Int((2*(1-f)*t^(-f)/(3*A*f*(1-1/(A*f*t^(f-1)))))^(1/2), t):
value(%) assuming A > 0, f > 0, f < 1;

(Medium-length result (1/3 of a screen) not worth transcribing here obtained with Maple 18.02.)

It isn't a bug. To aid in understanding this, you should remove the smooth option. This example, in a nutshell, shows the whole purpose of the range option to plots:-listdensityplot. Without the option, the range is 0..max(k, 1-k). For example, in a frame where k=0.75, the one cell with that value will be pure white, and the one cell with value .25 will have a grayscale value of .25/.75 = 33% of pure white. If you use range= 0..1, then the denominator is always 1, so .75 always appears as 75% of pure white and .25 always appears as 25% of pure white.

Two ways with unevaluation quotes:

seq('a[i], b[i]', i= 1..3);

''a[i], b[i]'' $ i= 1..3;

From your title, I assume that you mean that you want to plot a line with the geom3d package. But you also say "Line has coordinates x, y, and z." The geom3d:-line command has no way to specify the names of the coordinates, so I don't know what you mean.

Here's an example of how to plot in 3D a line from [1,1,1] to [2,2,2] using geom3d.

geom3d:-line(L1, [geom3d:-point(A,[1,1,1]), geom3d:-point(B,[2,2,2])]):
geom3d:-draw(L1, axes= boxed, orientation= [70,70], color= red);

What's your integration variable: t or x? The only way that you're going to get a numeric answer from this (which I assume is what you want---why else would you use evalf?) is if the integration variable is the same as the "free" variable in curva. If you want to maintain the generality of being able to specify the free variable at the time that the procedure is called, then you'll have to pass the free variable to both procedures. In other words, you passed t to curva, so you'll need to pass it to fun also. One way to do this that maintains most of the structure that you already have is to include t as the final return value of the sequence returned by curva and to make t the final parameter of fun. This will let you maintain the call structure fun(curva(m,t)). The code is

curva:= (m,t)-> (a-> (2*(cos(a), sin(a)*cos(a*m), sin(a)*sin(a*m)), t))(2*Pi*t);
fun:= (a,b,c,t)-> evalf(Int(add(_k*Pi*args[_k]*cos(_k*Pi*t), _k= 1..3)^2, t= 0..1));

And your example usage is

fun(curva(2,t));

     75.02998738

Update: The idea of passing around the free variable doesn't sit well with me. It seems artificial. Here's a formulation that doesn't need it.

curva:= proc(m)
local t, a:= 2*Pi*t;
     unapply~(2*[cos(a), sin(a)*~[cos,sin](a*m)[]], t)[]
end proc;

fun:= proc()
local k, t, a:= k*Pi;
     evalf(Int(sum(a*args[k](t)*cos(a*t), k= 1..nargs)^2, t= 0..1))
end proc;

fun(curva(2));

     75.02998738

It'd be better if you didn't include the output in your Questions until you learn how to properly format it. Just the input and error message are usually enough to get an Answer here on MaplePrimes.

Here's the relevant input and output upto and including the solve command.

restart:
v1 := int(cos(tau)*g(tau), tau = t0 .. t):
v2 := int(-sin(tau)*g(tau), tau = t0 .. t):
soln := C1*y1+C2*y2+v1*y1+v2*y2:
soln := combine(soln):
Eqs:= eval(soln, t = t0) = 0, eval(diff(soln, t), t = t0);
solve({Eqs}, {C1, C2});

So, those are your two equations (the second being considered to be equated to 0). The first has the variables C1 and C2. The second has neither. That pair of equations can't be solved for both C1 and C2. Notice that the solve command has given no output. That's because the equations can't be solved. So in your subsequent eval command, the % actually refers to the output of the last command which had non-NULL output, which is Eqs.

ex:= D[2,2](G)(x,y):
[op([0,0,..], ex)];

You can get the plot in your code if you reduce the precision of the numeric integration. This is done by changing int to Int and adding an epsilon option to the Int. Using epsilon= 1e-5 will make no visible difference in the plot. After you create the expression U, add the line

U1:= subsindets(U, specfunc(anything, int), J-> Int(op(J), epsilon= 1e-5));

Then simply use U1 instead of U in the plot command. You'll get this plot:

which shows that your iterative technique is far from the known exact solution.

In your second solve, one of the equations is redundant. With exact solutions, this wouldn't matter---solve would figure it out. But with floats, the redundant equation deviates enough from its exact value (for m > 3) that the system becomes inconsistent . I don't think that it's possible to solve this system exactly in any reasonable time (not sure about that). The solution is to seq the equations and variables for n= 3..m rather than n= 2..m.

There's another issue: Your second solve returns multiple solutions, even for m = 3. Your current code is ignoring all but the last of these solutions. In the code below, I return all solutions. Here's my version of your code with the two corrections. It works for m = 3, 4, 5. The m = 5 took about 20 minutes.

restart:
Digits:= 30:
m:= 3:  g:= 0.3:  nu:= 0.2:  a:= 1:
w:= unapply(add(b[n]*cos(n*r), n= 1..m), r):
W:= simplify(eval(w(r), solve((D@@2)(w)(1) + nu*D(w)(1), {b[1]}))):
d1:= diff(W,r):  d2:= diff(d1,r):
F:= int(((d2+d1/r)^2-2*(1-nu)*d2*d1/r)*r*(1+g*r/a)^3, r= 0..a)/int(d1^2*r, r =0..a):
Sol:= solve({seq(eval(diff(F,b[n]), r= n/m), n= 3..m)}, {b[n] $ n= 3..m});
FW:= map2(eval, F, [Sol]);

n:= 3:  N:= 5:
A:= Matrix((n,n), symbol= a):
B:= Matrix((n,n), symbol= b):
z:= Matrix(N-n,n):  Z:= Matrix(N,N-n):
<<A,z>|Z> + <Z|<z,B>>;