The reason for your print problems is a feature of the old-style matrices and arrays which isn't worth explaining here because you shouldn't be using them anyway. Your code can be vastly simplified and modernized to the new Matrices and Arrays like this:

restart:

Der:=proc(A,n)
local i,j,k,l,C;
     eval(
          Matrix((n,n), symbol= C),
          solve({
               seq(seq(seq(add(
                    A[i,j,k]*C[k,l] = A[k,j,l]*C[i,k]+A[i,k,l]*C[j,k],
                    k= 1..n), l= 1..n), j= i+1..n), i= 1..n-1)
          })
     )
end proc:

A1:= Array((1..2)$3, {(1,1,2)= 1}, storage= sparse):
Der(A1,2);

The above code for procedure Der is functionally equivalent to Tom Leslie's. If the Array A is highly sparse, as in the above example, there'll several occurences of the escaped local variable C. This may be a nuisance, and you may wish to replace it with a global variable. In that case, remove C from the local variables, and change the procedure's other four occurences of C to _C.

That is done by converting the first ODE's solution function for f(x) into something that'll return unevaluated when passed a symbolic argument and return numeric when passed a numeric argument and then declaring that as known in the dsolve command for the second ODE. Like this:

Sol_f:= dsolve({diff(f(x),x)=f(x), f(0)=1}, numeric):
F:= x-> `if`(x::complexcons, eval(f(':-x'), Sol_f(x)), 'procname'(x)):
Sol_g:= dsolve({diff(g(x),x)=F(x)*g(x), g(0)=1}, numeric, known= [F]);

Kitonum's Answer works also, but there are times when you want to solve the two as separate systems. And the above (using option known) will work for any numeric procedure, not just one that comes from a prior dsolve.

By the way, your title is perfect. I don't know what you see wrong with it.

That is surprising that it can't do it, because it can do it immediately if you replace the floats with symbols:

int(x*(1+a*x^4)^e*(C1-C2*x^4/(1+a*x^4)), x= 0..1);

     (1/8)*(4*C2*a^(5/4)*GAMMA(3/2+e)*GAMMA(1-e)*a^(-3/4+e)*hypergeom([1-e, -1/2-e], [1/2-e],                      -1/a)+4*C1*a^(1/4+e)*hypergeom([-e, -1/2-e], [1/2-e], -1/a) * a^(5/4)*GAMMA(3/2+e)*GAMMA(-e) * e-(4*        (1/2+e))*Pi^(3/2)*sec(Pi*e)*(C1*a*e+(1/2)*C2)) / (e*a^(3/2)*GAMMA(3/2+e)*(2*e+1)*GAMMA(-e))

evalf(eval(%, [C1= 2.906663106, a= 1.38, e= -0.4311594203, C2= 3.458929096]));

      1.00000000009749

This is the same (except for a little round-off error in the last digit) as what you'll get if you apply numeric integration (evalf(Int(...))) to your original integral.

If you want to compute an integral of a dsolve(..., numeric) solution, and one of the limits of integration is the initial point (or a boundary point) of the system, then the best way usually---the way that gives the best control of truncation error---is to build the integral directly into the system of ODEs, like this: Add to the system the ODE

diff(IntW(y), y) = W(y)

and the initial condition

IntW(ymin) = 0.

Then int(W(y), y= ymin..k) (which can't actually be directly computed like that) is equivalent to

eval(IntW(y), sol(k)); 

Edit: This idea comes from Joe Riel.

In the 2D input mode, when a variable is immediately followed by a left parenthesis, it is assumed that the variable is a function symbol being applied to the arguments inside the parentheses. In order to imply multiplication, you need to put a space between the variable and the parentheses. In other words, you need to change

p(4x^2+6) 

to

p (4x^2+6)

Oh, and with(student) serves no purpose in your worksheet. So why include it?

Use nested for loops, like this:

S:= 0:
for m from 1 to 5 do
     for j from 0 to m-1 do
          k:= m-1-j;
          S:= S + 2^m*(D@@j)(f)*(D@@k)(f)
     end do
end do;

where f(j,k) is the summand.

Or, do you really need k? This works:

add(add(2^m*(D@@j)(f)*(D@@(m-1-j))(f), j= 0..m-1), m= 1..5);

Yes, you are essentially right about the problem being due to extra brackets. The cause of the brackets is that XACT and YGPA are created as two-dimensional Arrays with a degenerate second dimension. Essentially, they're 120 x 1 matrices whereas you're probably thinking of them as vectors. XACT[9], for example, extracts the ninth row, not the ninth element. Those are essentially the same thing---the rows only have one element each---but to Maple your L1 looks like a list of pairs of 1x1 matrices rather than a list of pairs of numbers. We could convert XACT and YGPA to Vectors, and if I was going to work with them at length, I'd probably do that. But for this job that's not necessary. And they don't need to be converted to lists either. They can be plotted as they are, like this:

F:= Statistics:-Fit(a+b*x, XACT, YGPA, x);

Scatter:= plot(<XACT|YGPA>, style= point, symbol= circle, color= blue, symbolsize= 9):
Line:= plot(F, x= (min..max)(XACT), thickness= 2):
plots:-display(
     [Scatter, Line],
     title= "ACT vs GPA",
     labels= ["ACT", "GPA"],
     ytickmarks= [1,2,3,4],
     view= [14..35, 0.5..4.0]
);

The command to combine plots is plots:-display.

By the way, just out of curiosity, are those GPAs measured at high-school graduation, college graduation, or some other time?

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;