Why not just divide by y?

expr/~y;

Another way, which'll work for any number of lists of the same size, is f~(A, B).

There's no way to do it using map. (I was wrong: Kitonum found a way with map, shown below.) There is zip, as shown by Acer, but it won't work for more than two lists.

If you don't like ~, you could define a procedure so that you only need to use ~ once:

MapThread:= f-> f~(_rest):

Then MapThread(f, A, B) will work for all future uses.

Before you can understand how that example works, you need to understand how divide-and-conquer recursion works, without the use of multiprocessing. The procedure below multiplies a list of numbers. I chose this instead of adding a list of numbers because the below turns out to be a very efficient method for multiplying. The algorithm logic would be identical for adding. So, please study this procedure closely, and let me know if you have any questions about it.

The keyword thisproc invokes a recursive call to the procedure.

Mul:= proc(L::list, i::posint:= 1, j::posint:= -1)
local k;
    if j = -1 then thisproc(L, 1, nops(L)) #initial call
    elif i = j then L[i]
    else
        k:= iquo(i+j, 2); #floor((i+j)/2)
        thisproc(L, i, k) * thisproc(L, k+1, j)
    fi
end proc
:

The above is a very efficient procedure for multiplying integers. It beats the default builtin mul by a huge factor:

n:= 2^17-1:  L:= ['rand'()$n]:
r1:= CodeTools:-Usage(Mul(L)):
memory used=107.55MiB, alloc change=0 bytes, 
cpu time=1.78s, real time=1.64s, gc time=484.38ms

r2:= CodeTools:-Usage(mul(L)):
memory used=38.42GiB, alloc change=0 bytes, 
cpu time=75.06s, real time=51.87s, gc time=64.33s

r1 - r2;
                               0

If the problem is as you describe, then a workaround is to put a page break (Ctrl-Enter) after the plot.

Click on the <> icon on the toolbar. Paste or type text in the white dialog box that appears. When you dismiss the dialog, the grey box will appear, and you can continue to edit it. The text will not appear correctly formatted until it's in the grey box.

Try just closing worksheet A before running B.

Procedure parameters can't have indices, such as beta[1]. They can have subscripts, as in beta__1.

The legal values for symbol are "asterisk", "box", "cross", "circle", "diagonalcross", "diamond", "point", "solidcircle", "solidbox", or "soliddiamond", with or without quotes.

The best way to get a better curve is to have more data points. Is that possible?

The situation is not as described by Rouben. A complete symbolic solution is mathematically trivial (although too big to do by hand). I don't know why dsolve won't do it. Here's how to do it. After your A1 is entered as above, do

Eq7:= value(dsolve({A1, f[3](0)=0, D(f[3])(0)=0})):
BC1:= eval(rhs(Eq7), x= 1) = 0:
BC2:= eval(diff(rhs(Eq7), x), x= 1) = 0:
Csol:= solve({BC1, BC2}, indets(Eq7, suffixed(_C))):
Eq7A:= eval(Eq7, Csol):

The solution Eq7A is easily displayable and about half the size of A1.

What Rouben failed to see was that your ODE was of the trivial form diff(f(x), x$4) = P(x), where P(x) is just a polynomial (with coefficients with 5 parameters).

Your use of seq is not bad. All of the empty lists can be removed at once by subs([]= (), ...). Here's my procedure:

Tlist:= (L::listlist)-> 
    subs([]=(), [seq]([seq](x[i], x= L), i= 1..nops(L[1])))
;

This matrix-transpose-based version might be faster. I haven't speed tested it:

Tlist2:= (L::listlist)->
local _; 
    (M-> [seq]([seq](M[i]), i= 1..nops(L[1])))
        (subs(_= (), rtable(subs([]= _, L)))^+)
;

If you want to apply an operation to the terms of an expression, there's no need to use lists or convert.  All that you need is

MapTerms:= (f,e)-> `if`(e::`+`, map(f,e), f(e)):

Example:

MapTerms(f, a+b+c);

The number of boundary conditions required is the differential order plus the number of unspecified parameters. It often happens that a symbol is mistyped and it gets interpretted as a parameter.

Regarding your followup question: Your Eq3 can be easily solved for omega and then plotted with the regular plot command:

plot([solve(Eq3, omega)], R);

Standard coeff works:

coeff(tt, N__s, -3)

The problem is that you had 1/N__s rather than N__s as the second argument.

See the fnormal command (help page ?fnormal). Note that it's necessary to follow it with simplify(..., zero) to remove the Is.