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.

The problem is that when you load a Units subpackage (such as with(Units:-Simple)), all of the basic arithmetic operations are replaced with units-aware versions. This is one of the hazards of the with command. The replaced operations are listed as part of the output of the with command. One of those operations is max. I don't know why the Units version doesn't work in your case, but you can call the original max with :-max. So, you can use :-max~(A, B), where A and B are matrices, or you can avoid the need to create a zero Matrix by using :-max~(0, B). The :- prefix can be used on any name to guarantee that it refers to a global version.

Perchance, have you created your own units? I didn't look at the code underlying your embedded components because I don't know how to.

A homogenous linear constant-coefficient recurrence can be done quickly by matrix powering:

((<0,1,0; 0,0,1; 1,-2,1>^(50-3)).<2,-1,0>)[-1];

The last row of the matrix is the coefficients. The other rows are always a staircase pattern of elementary unit vectors as shown. The 3 is the number of coefficients. The vector is the initial conditions. The -1 selects the last entry of the result.

If the exponent is large, say 100,000, then this method is significantly faster than methods that require computing the entire preceding sequence of terms.

A:= proc(n::And(posint, Not(1)))
local a, a1:= 1.75, a2:= 2.25;
    to n-1 do 
        a:= sqrt(a1+4)/6/a2;
        (a1, a2):= (a, a1)
    od;
    a
end proc
:
A(182);

This avoids the need to store the entire sequence of previously computed values.

Use a procedure, like this:

MySum:= (j::integer)->
    add(combinat:-numbcomb(j+i, i+1)*w^(j-1)*(1-w)^(i+1), i= -1..j-2)
: