The command that you want is coeff rather than coeffs:

coeff(r, x, 1);

will return the coefficient of the first power of x. You may then pass this to solve.

Here's the much faster way that I promised---the way that doesn't need to trace back over the path to determine that it'll self intersect on the next step. This generates 1024 walks in an eighth of a second on my machine, the memory usage is insignificant, and the code is much simpler. Indeed, there are no if statements and a only single for loop of two statements.

(Upon re-reading your Question, I see that you're using the word "line" to refer to what the rest of the English-speaking mathematical world (AFAIK) refers to as a random walk. To avoid confusion, I use "walk" in the following. You may read this to yourself as "line". The term path is reserved (by the English-speaking mathematical world) for walks that don't intersect themselves.)

restart:

(M,N):= (2^10, 2^7): #(num walks, max length of walk)

#Tom Leslie's "turn" functions, simplified slightly:
gostraight:= (p,q)-> 2*~p-~q:
turnLeft:=     (p,q)-> p+~`if`(p[1]=q[1], [q[2]-p[2], 0], [0, p[1]-q[1]]):
turnRight:=   (p,q)-> p+~`if`(p[1]=q[1], [p[2]-q[2], 0], [0, q[1]-p[1]]):
dirs:= [gostraight, turnLeft, turnRight]:
R3:= rand(1..3):

OneWalk:= proc(MaxLength::posint)
local
     W:= table([1= [0,0], 2= [1,0]]), #The walk
     V:= table([[0,0]= NULL]), #The points visited, without the terminus
     k
;
     for k from 3 to MaxLength+1 while not assigned(V[W[k-1]]) do
          V[W[k-1]]:= NULL; #Mark previous point as visited.
          W[k]:= dirs[R3()](W[k-1], W[k-2]) #Generate next point.
     end do;
     convert(W, list)
end proc:

Walks:= CodeTools:-Usage(['OneWalk(N)' $ M]):

memory used=17.46MiB, alloc change=32.00MiB, cpu time=125.00ms, real time=125.00ms, gc time=0ns

They can be plotted by simply plot(Walks); however, 1024 is probably too many to appreciate on a plot, so you can do, for example, plot(Walks[100..150]);

I've found all the roots in one period of the function in under a second. Since the function is periodic, this is essentially the same as finding all the real roots.

First, I do a little simplification with the combine command.

restart:
Digits:= 15:
f:= 2*cos(0.5*x)*sin(0.5*x)*cos(3.775*x) +
    2.2075*cos(0.5*x)^2*sin(3.775*x) -
    0.453*sin(0.5*x)^2*sin(3.775*x)
:
fc:= convert(combine(f), rational);

Several things are obvious from this simplified form (several of which were also obvious from the unsimplified form):

1. The function is odd (in the symmetry sense), so the negative of every positive root will also be a root.
2. Zero is a root.
3. The function is periodic with period 80*Pi.
4. The function is well-behaved. It is entire (no discontinuities) and finding roots should be easy.
5. Items 1 and 3 imply that finding the roots on 0..40*Pi will be sufficient.

My first choice for finding all the real roots of such a function is RootFinding:-NextZero.

ffc:= unapply(fc,x):
R:= table(): #To store the roots
r:= 0: #First root
st:= time():
while is(r < 40*Pi) do
     R[r]:= [][]; #Right side can be anything.
     r:= RootFinding:-NextZero(ffc, r, guardDigits= 3);
     if not r::numeric then
          print(r);
          break
     end if
end do:
time()-st;

     0.670

R:= {indices(R, 'nolist')}; #Make table into a set.

(... large set of roots omitted...)

nops(R); #Count the roots.

     191

max((abs@ffc)~(R)); #Check residuals

Since I used three guardDigits, these roots are only guaranteed accurate to 12 digits.

R:= evalf[12](R):

There's no single command to do a double break, but it's a fairly common operation. Here's a standard way of handling it:

for i to M do
     broke:= false;
     for j to 40 while not broke do
          R[i,j];
          for k to j-2 do
               if R[i,j]=K[i,k] then
                    x:= j;
                    print(b[i]=j);
                    broke:= true;
                    break
               end if;
               L[i,j]:= R[i,j]
          end do;
          #The following line isn't needed unless there's code in the j
          #loop after the end of the k loop.
          if broke then break end if;
          # ...possible additional j-loop code.
     end do
end do;

For aesthetic reasons, I usually think of some way to avoid the double break, but I can't state a standard way to avoid it in all cases. If aesthetics are of no concern, then I see no reason to not do it. (But see my "Better ways" Reply below.)

Simply,

Gcd(f,g) mod 3;

Note that the capital G distiguishes this from the "regular" command gcd. See the help at ?mod. This will have a list of all the commands designed to work with mod over finite rings.

What's immediately obvious to me is that you're not updating X and Y. I think that you think that your updates of R implicitly also update X and Y, but that doesn't happen. They can only be updated when they appear on the left side of an assignment statement.

My advice is to not use R at all; just maintain X and Y. That is, for each of the three values of r there should be two assignment statements: one for X and one for Y.

This is one of the most common Maple errors, called premature evaluation. The functions I1 and I2 obviously can't be evaluated numerically until t has a numeric value. In the plot command you invoke I1(t) and I2(t). These invocations are evaluated before t is given any numeric value. One solution is this:

plot([I1,I2], 0..200, color= [red, blue]);

P.S.: The above is a correction to the premature evaluation problem only. I haven't tested whether this particular plot actually works. Markiyan got a nonreal value. If you get too many nonreal values, the plot command will give you an empty plot.

Here is an example random problem with 15 variables, nine general constraints, and all variables restricted to the unit interval:

V:= [x||(1..15)]:
Objective:= randpoly(V, degree= 1, dense):
Constraints:= ['randpoly(V, degree= 1) >= 0' $ 9]:
Optimization:-LPSolve(Objective, Constraints, (V =~ 0..1)[]);

This code executes instantaneously. However, you specified the open interval 0 < x < 1. Is that what you really want, as opposed to the closed interval 0 <= x <= 1? Open-interval constraints can't be done by the program, and it isn't even clear whether they can be solved at all. Continuous functions are guaranteed to have extrema on closed intervals, not on open intervals. You can get an approximate solution by choosing a small positive number eps and using the closed interval eps..1-eps. I don't know off the top of my head whether the limit of the solutions thus obtained will necessarily converge to the true solution as eps approaches 0; my guess is yes, but that there's no pratical way to compute how small eps needs to be.

Your procedures F and L are correct, although L's performance can be improved significantly by including option remember.

The procedure S is very wrong, and it's not entirely clear how to correct it. I think that you may mean this:

S:= proc(n::nonnegint)
local i;  
     for i from 1 by 1 to n do
          printf("%a\n", L(i)^(2)-5*F(i)^(2))
     end do  
end proc;

1. There's no need to refer to S inside S.
2. You were missing end do.
3. Formatted printing is done with printf, not print.
4. When using printf, you need to handle all spacing and line breaks. The \n is a "new line".

Your L can be shortened to

L:= proc(n::nonnegint)
option remember;
     if n = 0 then 2 elif n = 1 then 1 else L(n-1)+L(n-2) end if
end proc;

or

L:= proc(n::nonnegint)
option remember;
     if n = 0 then 2 elif n = 1 then 1 else thisproc(n-1)+thisproc(n-2) end if
end proc;

When you refer to a procedure from within it's own body, you can refer to it as thisproc. This is a tiny bit faster than referring to it by name, and it also means that you won't have to change the procedure's body at all if you decide to change the procedure's name.

Example:

P1:= plot3d(x*y, x= -1..1, y= -1..1):
P2:= plots:-pointplot3d([[0,0,0]], symbol= solidsphere, symbolsize= 32, color= red):
plots:-display([P1,P2]);

Why is your data organized in this strange way? Something tells me that you're making it more complicated than it needs to be. If each letter is to be paired with a number, then store them as pairs:

map(L-> {seq(L[2*k-1]=L[2*k], k= 1..nops(L)/2)}, data);

You see that when you use this more natural format, the sets are automatically sorted by letter. Unlike the other two Answers, mine will work on sublists of any length.

Your manual solution is faulty because the Laplace parameter s still appears in the equation. Since no gamma appears in your solution, I will assume that the coefficient gamma in the second original equation was supposed to be lambda.

These can be solved easily with Maple, like this:

restart:
Ps:= [P[0],P[1]]:
eqs:= Ps(s) =~ [1/(s+lambda)+mu*P[1](s)/(s+lambda), lambda*P[0](s)/(s+mu)]:
Ls:= solve(eqs, Ps(s))[];

Ps(t) =~ inttrans[invlaplace]~(rhs~(Ls), s, t);

If the step size is (b-a)/n, then including a and b, there are n+1, not n, evenly spaced points in the interval.

They can be generated like this (I take the name linspace from the equivalent Matlab function):

linspace:= (a,b,n)-> [seq(a+(b-a)*k/(n-1), k= 0..n-1)]:
linspace(0,1,50);

The procedure f can be applied to them like this:

f~(linspace(0,1,50));

While time is useful for testing a block of code (which should be in a single execution group), it is better to use CodeTools:-Usage to test a single procedure call. It returns more-detailed information. Example:

P:= `*`('randpoly(x, degree= 2^9, dense)' $ 2^7):
Prod:= CodeTools:-Usage(degree(expand(P)));
memory used=0.58GiB, alloc change=257.00MiB, cpu time=16.59s, real time=16.62s, gc time=0ns

     Prod:= 65535

By default, Usage returns directly the result of the computation and returns its usage report as a side effect. This default can be changed with Usage's output option.

This is a common source of bugs in Maple code: The code is expecting an expression sequence. If the sequence has only one member, then indexing the sequence becomes indexing of the member, which is not what was intended. The solution is to turn the expression sequence into a list by enclosing it in square brackets:

Answers:= [solve(equation that needs solving)];

This will work even if solve returns nothing. Now if Answers has only one solution, Answers[1] will return that solution in its entirety.

Terminology: Neither a list nor a sequence should be called a matrix, even when it can be double-indexed. And an expression sequence shouldn't be called a list, as in the title to this Question.

The number of elements of a list is returned by nops. This'll tell you what are the safe indices to use.

You should not use indexing to extract the a and b values from an individual solution. Use eval instead, as in

eval(a, Answers[1]), eval(b, Answers[1]);