One problem is your method for specifying the boundary conditions. The first should be specified as

 Phi(1)=1

The second is problematic. Are you trying to specify a boundary condition at a symbolic value, p1? I don't believe that numeric dsolve will handle that (but am not sure). Another problem is the existence of the complex values in the boundary condition. There are also problems with some of your assignments. For example

V:=r->V0(r)+f(r)^2*mu^2-(omega-m*Eta)^2+(f(r)^2/(r^2))*(l(l+2*N)-m^2*(1-(r^2/(h(r)^2)))+4*(1-sigma)*(h(r)^2/(r^2)-1));

is suspect, the term l(l+2*N) is probably not what you intend (it always evaluates to 2, since l is assigned 2:

 l := 2: l(l+23+a);
                                       2;

Probably you meant l*(l+2*N). The error could be just the cut/paste from the standard gui, but it is usually best to use explicit multiplication.

The problems are the references to A[k,j] and M[k,j]. You cannot reference an rtable (a data structure that encompasses Matrices and Vectors) with a symbolic index. For example:

V := Vector([1,2,3]):
sum(V[k], k=1..3);
Error, bad index into Vector

If the sum is changed to an add, then the problem goes away because an integer value is substituted for the k before it is evaluated (add works at a lower level than sum):

add(V[k], k=1..3);
                                       6

In your example, changing sub to add is useful, regardless, because add is usually better for summing over a given range. However, that alone won't fix the problem because the j term is symbolic. To fix that just wrap the whole expression in your do loop:

for k from 1 to 5 do
    for i from 1 to 10 do
        u[i+1][k] := u[i][k]-w*(1/(M[k,k]))*((add(M[k,j]*u[i][j],j=1..5))
                                             +e*b[k]+(add(A[k,j]*p[j],j=1..5))
                                             -(add(M[k,j]*u[i+1][j],j=1..k-1))
                                             +(add(M[k,j]*u[i][j],j=1..k-1)));
    end do;
end do:

Now Maple is happy. Alas, that doesn't fix all the problems. The handling of u is mostly likely wrong (I'm not sure what you intend).

LegendreP(1,1,cos(theta));
                         LegendreP(1, 1, cos(theta))
simplify(%);
                                 (1/2)                 (1/2)
                 (cos(theta) - 1)      (cos(theta) + 1)     
expand(%^2);
                                             2
                              -1 + cos(theta) 
simplify(sqrt(%),symbolic);
                                I sin(theta)

Note that you should change the relations to

0 ≤ x ≤ a, a < x ≤ b

Otherwise there is an indeterminancy when x = a (which Maple resolves, but is poor form).

Much as I enjoy doing homework, I prefer doing my own. However, I'll give you a few hints. First, I suspect that the asterisk used in [g*f] is intended to be the composition operator, not multiplication. Use @ (see ?@), the Maple composition operator, to do this, after assigning f and g as functional operators (see ?operators,functional). Here is an example of composing a function with itself

f := x -> x^2 + 2:
h := f@f:
h(x);
           (x^2+2)^2 + 2

Note that one can also do

h := f*f:
h(x);
            (x^2+2)^2

so the proper answer depends on the interpretation of the asterisk.

Use plotsetup:

plotsetup('gif', 'plotoutput' = "/tmp/my_plot.gif"):
plot(sin, 0..Pi);

To do this in a loop, you might try

for i to 5 do
   plotsetup('gif', 'plotoutput' = sprintf("/tmp/my_plot%d.gif", i));
   plot( f[i], ... );
end do:

where f[i] is the expression to be plotted.

Here's the comparison without the algorithmic optimizations.

               +-------- timing -------+--------- memory ---------+
Procedure      | maple (s) | total (s) | used (MiB) | alloc (MiB) | gctimes
---------------+-----------+-----------+------------+-------------+--------
tp_original    |    22.69  |    23.08  |   246.61   |     36.43   |    6
tp_evalhf      |     2.06  |     2.16  |     0.00   |      0.00   |    0
tp_autocompile |     0.15  |     0.38  |     0.55   |      0.00   |    0
tp_compile     |     0.09  |     0.09  |     0.00   |      0.00   |    0
---------------------------------------------------------------------------

tp_compile looks like

tp_compile := Compile:-Compile( proc() ... end proc):

It doesn't run any faster than tp_autocompile, the timing difference is that the compilation time is not included in the run. To measure those accurately I'd need to execute them several times. Declaring the local counters as integers doesn't have any measurable effect (measured with multiple executions).

You can also use the autocompile option, which is frequently faste than calling evalhf. Also, there are two small optimizations you can make, not really pertinent to the broader topic. First, don't compute the square root, just compare against 100.0^2. Second, break out of the inner loop if Temp3 ≥ 100.0^2. Here are the revisions:

testproc3 := proc()
local tp;
    tp := proc()
    local i, j, k, total, Temp3, Temp1, Temp2;
        total := 0;
        for i to 100 do
            Temp1 := i*i;
            for j to 100 do
                Temp2 := Temp1+j*j;
                for k to 100 do
                    Temp3 := Temp2+k*k;
                    if Temp3 >= 10000.0 then
                        break;
                    end if;
                    total := total+1;
                end do
            end do
        end do;
        total/0.1e6
    end proc:
    evalhf(tp());
end proc:


testproc4 := proc()
local i, j, k, total, Temp3, Temp1, Temp2;
option autocompile;

    total := 0;
    for i to 100 do
        Temp1 := i*i;
        for j to 100 do
            Temp2 := Temp1+j*j;
            for k to 100 do
                Temp3 := Temp2+k*k;
                if Temp3 >= 10000.0 then
                    break;
                end if;
                total := total+1;
            end do
        end do
    end do;
    total/0.1e6
end proc:

With testproc1 and testproc2 being yours and Jacque's procedures, a comparison of the four gives

            +-------- timing -------+--------- memory ---------+
Procedure   | maple (s) | total (s) | used (MiB) | alloc (MiB) | gctimes
------------+-----------+-----------+------------+-------------+--------
testproc1   |    22.97  |    24.74  |   247.12   |     37.68   |    6
testproc2   |     2.06  |     2.05  |     0.00   |      0.00   |    0
testproc3   |     1.04  |     1.03  |     0.00   |      0.00   |    0
testproc4   |     0.27  |     0.55  |     3.27   |      2.56   |    0
------------------------------------------------------------------------

Slightly simpler is just

map(sin, [3,4,5]);

or

evalf(map(sin,[3,4,5]));

Another way to do this is

[seq(sin(x), x=[3,4,5])];

It builds a sequence and encloses it in a list.

The rest of your post is missing. If it includes a less-than sign, substitute the html equivalent `<' (not including the quotes). Note that the range in the seq assignment to datalist should be 1..n, not 0..n. Also where is plot2?

The way to do this is to compute the compute the intercept point (use yint := (ty*x-y*tx)/(x-tx)) then plot/write that value. I'm not sure what you want to plot it against, maybe just the index? Better would probably be against the distance of the incoming ray from the central axis. Here's how to do that (I added the x_vs_yint table to store the values; using a table, while not necessary here, is more efficient than building a list a term at a time):

with(plots):
with(plottools):
mirror := x -> 8 - sqrt(64 - x^2):     # Mirror
slope := x -> x/sqrt(64 - x^2):        # Slope at each point on mirror
width := 5.0:                          # x = -slope..slope
height := 2 * width:                   # y = 0..height
n := 30:                               # number of light rays

unassign('x'):
unassign('y'):

blk := plot(mirror(x),
            x = -width..width, y = 0..height,
            axes = NONE, scaling = CONSTRAINED, color = black):

rays := {}:
x_vs_yint := table();
cnt := 0;

for i from 0 to n do
    x := -width + 2 * i * width/n:
    y := mirror(x):
    rays := rays union {line([x, height], [x, y])}:   # Incoming ray
    dr := slope(x):
    if (abs(dr) > 0.001) then          # if slope = 0 no outgoing ray
        m := (dr^2 - 1)/(2 * dr):
        ix := (height - y)/m:
        if ix < -width then
            tx := -width:
            ty := y + m * (tx - x):
        elif
        ix > width then
            tx := width:
            ty := y + m * (tx - x):
        else
            ty := height:
            tx := (ty - y)/m + x:
        fi:
        yint := (ty*x-y*tx)/(x-tx);
        cnt := cnt+1;
        x_vs_yint[cnt] := [x,yint];
        rays := rays union {line([x,y], [tx, ty])}:  # outgoing ray
    fi:
od:
rd  := display(rays, color = red, axes = NONE, scaling = CONSTRAINED):
unassign('x'):
unassign('y'):
display({blk, rd});
plot([seq(x_vs_yint[i], i=1..cnt)]);

Do not use with inside a procedure. It is for top-level use only. To access components of a module from a procedure or module, use the use ... end use statement or the uses phrase. You can also give the full path to the command. For example

myproc := proc(s::string)
uses StringTools;
   IsAlpha(s);
end proc:

By default, LinearAlgebra uses, when possible, hardware floats to do Matrix computation of floats. You can explicitly turn this off by setting the environmental variable UseHardwareFloats false.

M := <<0.5|0.5>,<0.3|0.7>>;
                                    [0.5    0.5]
                               M := [          ]
                                    [0.3    0.7]
M^2;
                [0.400000000000000024    0.599999999999999978]
                [                                            ]
                [0.359999999999999986    0.639999999999999904]
UseHardwareFloats := false:
M^2;
                               [0.400    0.600]
                               [              ]
                               [0.360    0.640]

Another possibility is to use the older linalg package, but you need to convert to matrices.

m := convert(M,matrix):
evalm(m^2);
                                [0.40    0.60]
                                [            ]
                                [0.36    0.64]

How about posting a snippet of code that illustrates the problem. You also might want to try using Array rather than array, but its hard to say without seeing what you are doing.

A better approach would be to delete them all in one call.

M := DeleteRow(M, [seq(`if`(M[i,2]=M[i,3],i,NULL), i=1..RowDimension(M))]);

Another way to do this is to keep the rows that are "good".

M := M[[seq(`if`(M[i,2]=M[i,3],NULL,i), i=1..RowDimension(M))], [1..-1]];

To insert a less-than sign in html use <.