In your second worksheet there is a procedure FiveSpeedGearBox_R which contains the two lines

eq:= piecewise(ig=1,i[1],ig=2, i[2],ig=3,i[3],ig=4,i[4],ig=5,i[5],1);
return eq(ig);

The expression eq(ig) is gibberish. Unfortunately, Maple does not recognize it as such and lets you run with it. A piecewise expression is not a procedure---it does not take an argument. Replace these two lines with the single line

return piecewise(ig=1,i[1],ig=2, i[2],ig=3,i[3],ig=4,i[4],ig=5,i[5],1);

After making that change, the solve runs for about 90 minutes and returns a solution without error or warning. That solution is five screens "wide" by 56 screens "tall". Here is the executed worksheet:

withGB.mw

You can use a recursive procedure.

X:= proc(Xo,Ro,n)
option remember;
     if not n::nonnegint then return 'procname'(args) end if;
     if n=0 then return Xo end if;
     Ro*thisproc(Xo,Ro,n-1)/(1+thisproc(Xo,Ro,n-1))^4
end proc:

BTW, I corrected your spelling of Hassell.

You must use N:= floor(a/h) or trunc(a/h). If you use round(a/h) then there's a risk that your x's will get larger than a. Note that trunc is semantically equivalent to floor for nonnegative numeric arguments, but trunc is simpler and faster.

test:=module()
uses DT = DocumentTools;
export abcd,tt,vv;
abcd:=proc()
    tt:=2000;
    vv:=DT:-GetProperty('Slider0', 'value');
    DT:-SetProperty('Label0', 'caption', vv);
end proc;

end module;

Your module test, which I've copied above, contains a procedure abcd which sets the values of some of the exports of the module. However, abcd is merely defined; it is never executed. So, right above the end module put the line abcd().

LinearAlgebra:-GramSchmidt({< 22,11,5 >, < 13,6,3 >, < -5,-2,1 >});

The imaginary parts are due to round-off error. Note that they are quite small. Their exact value is zero. It is impossible to avoid these imaginary parts if you solve a cubic by the cubic formula and and then substitute floating-point values for the coefficients. A plot shows that the cubic polynomial has three real roots. So, just throw away the imaginary parts. You can apply the Re command to each solution.

You want to use an Array of Records. See ?record .

I suggest that you enter your graph the normal way, by using the command GraphTheory:-Graph. GRAPHLN is intended for internal use. If you entered it the normal way, then the error pointed out by Kitonum would not have happened.

restart:
Eckert:= 0.5:
sigma:= 15:

eq1:= diff(f(eta), eta)^2-f(eta)*diff(f(eta), eta$2) =
diff(f(eta), eta$3)+k_1*(2*diff(f(eta), eta)*diff(f(eta), eta$3)
+diff(f(eta), eta$2)^2-f(eta)*diff(f(eta), eta$4)):

N:= 10:
bcs1:= f(0) = 0, D(f)(0) = 1, D(f)(N) = 0, (D@@2)(f)(N) = 0:

eq2:= diff(theta(eta), eta$2)+sigma*f(eta)*diff(theta(eta), eta) +
sigma*Eckert*(diff(f(eta), eta$2)^2+k_1*diff(f(eta), eta) *
diff(f(eta), eta$2)^2-f(eta)*diff(f(eta), eta$2)*diff(f(eta), eta$3)) = 0:

bcs2 := theta(0) = 1, (D(theta))(0) = -2.876:

Sol:= proc(k1)
     if not k1::realcons then return 'procname'(args) end if;          
     dsolve(
          eval({eq||(1..2), bcs||(1..2)}, k_1= k1),
          numeric, method= bvp[middefer], maxmesh= 2^13, abserr= 1e-4
     )
end proc:

plots:-animate(
     plots:-odeplot,
     [Sol(k1), [[eta, f(eta)], [eta, theta(eta)]], eta= 0..2],
     k1= 1.5..2.5
);

Yes, you can know the order of the error. The order of the error is an inherent property of the method; it does not depend on whether the exact solution is known.

Using solve instead of fsolve will give both solutions in this case, although I can't say that this is generally true. Certainly fsolve always gives at most one solution to a transcendental equation.

n:= LinearAlgebra:-RowDimension(A):
P:= (LinearAlgebra:-IdentityMatrix(n)-A)^(-1).A;

The syntax requires

dsolve({bs} union sys1);

Note that you used {} on sys1 and not on bs. You need to take that into account in your call to dsolve.

But there is a deeper problem. The system cannot be solved for the last constant. You can impose the first five boundary conditions like this:

dsolve({bs[1..5]} union sys1);

This will give you a system with one constant. For now, I leave it to you or someone else to address why it cannot be solved for this constant using your conditions.

This will be easiest to do in the case where the curves in the original are all one color and you want to convert them to all one color. Let's suppose that you want to convert them all to black. Then do

macro(C= ColorTools):
subs(
     indets(P, specfunc(anything, COLOUR))[] =
     COLOUR(RGB, C:-RGB24ToRGB(C:-NameToRGB24("black"))[]),
     P
);

Note carefully that COLOUR is spelled with a U and that ColorTools is spelled without a u (a result of a Canadian identity crisis :-) ).

If you are dealing with a case of multiple colors, something can still be done---let me know.

Use

return fff, ggg;

or

return [fff, ggg];

Use the former if you are prepared to accept multiple outputs. Use the latter if you want to treat the multiple outputs as a single entity.