To turn m into a matrix you could do

M:=Matrix(op~([entries(m)]));

There are several strange errors, e.g. in the following

while t do

if t >= 0 and t if t > 9.5 and t if t > 12.875 and t if t > 15.5 and t if t > 20.75 and t

What did you mean?

Another thing:

p.t := display({p_0 + V_1 / V_0 * K})

doesn't make any sense. Plots can be displayed, expressions cannot.

I think the problem lies in what you don't show us. You don't tell us what phi, a, b, f, or epsilon[step] are.

And if or when you do, please don't give us a picture, but text so we can try the code ourselves. Instead of text you could upload a worksheet.

If the order of the ode is n and you want a conventional solution which is n times differentiable on the real line you only have the n+1 equations f^(k) (0)=g^(k) (0) where k = 0, 1, ..., n.

But maybe you could give us the code as text, alternatively upload a worksheet.

In procedures evaluation takes place only once. So in deq  X and Y are not evaluated to what you want. An easy remedy is to move the definitions of X and Y up before the definition of deq.

But what are XX and YY? Should they just be X and Y?

Finally, you ought to include the line

uses plots;

right after the declaration of locals since you are using display from the plots package inside the procedure.

Since you have made the assignment zeta := tau-c*z; you must use another name.

Using another name (here zeta1) and still with zeta assigned to tau-c*z you can do as follows,

T1:=subs(zeta=zeta1,T);

dsolve(T1);

As an alternative, start out like this where you avoid assigning to zeta:
U := V(tau-c*z)*exp(I*(k*tau-w*z));
Then at the end you can do

T1:=subs(tau-c*z=zeta,T);
dsolve(T1);

Here are two ways. The first uses FindMinimalElement in the ListTools package, the other a double loop.
In both cases R contains the row numbers of the elements.

N:=10:
A:=LinearAlgebra:-RandomMatrix(N,generator=0..200);
R:=[seq(ListTools:-FindMinimalElement(convert(abs~(A-~100)[..,c],list),position)[2],c=1..N)];
A[R[3],3]; #In the third column the element in row number R[3] is closest to 100.

## The second method
R:=Vector(N):
for c to N do
   rm:=1;
   for r from 2 to N do
       if abs(A[r,c]-100)<abs(A[rm,c]-100) then rm:=r end if
   end do;
   R[c]:=rm
end do:
R;
A[R[3],3];

If the matrix is e.g. 100x100 or even larger you will notice that the second method is a lot faster.

You most likely don't get a more explicit symbolic result than what you got since it would involve solving a polynomial of degree higher than 4.

However, the result is not useless as can be seen from the example below.

(Assumptions don't help any).

restart;
P:=(b8*Omega^8+b6*Omega^6)/(a10*Omega^10+a8*Omega^8+a6*Omega^6+a4*Omega^4+a2*Omega^2+a0);
Omega_r1:=beta*Omega_c;
Omega_r2:=(2-beta)*Omega_c;
P_I:=int(P,Omega=Omega_r1..Omega_r2);
#Example
param1:=map(rcurry(`=`,1),indets(P,name) minus {Omega});
param2:={beta=1/2,Omega_c=1};
eval(P_I,param1 union param2);
evalf(%);
simplify(%);

#Compare this to a numerical integration (notice the capital I in 'Int'):
evalf(eval(Int(P,Omega=Omega_r1..Omega_r2),param1 union param2));

Use 'dsolve/numeric', 'animate', and 'pointplot'

restart;
with(plots):
g11:=1;g21:=g11;

z:=exp(I*((t/(2*Pi))+(Pi*j)));

f1:=evalc(Re((I/(2*Pi))*((g11/(subs(j=1,z)-(x(t)+I*y(t))))+(g21/(subs(j=2,z)-(x(t)+I*y(t)))))));

q1:=evalc(-1*Im((I/(2*Pi))*((g11/(subs(j=1,z)-(x(t)+I*y(t))))+(g21/(subs(j=2,z)-(x(t)+I*y(t)))))));

sys1:={diff(x(t),t)=f1,diff(y(t),t)=q1}; #Using a set

res:=dsolve({op(sys1),x(0) =0, y(0) = 1},numeric,output=listprocedure);
X,Y:=op(subs(res,[x(t),y(t)]));
animate(pointplot,[[[X(t),Y(t)]],symbol=solidcircle,symbolsize=15], t=0..500,frames=1000);

Maybe you could use a try ... end try statement:

ReInit := proc(dsol,tf::realcons)
    local imax,eps,i, pts, ini ;
    imax:=100;
    eps:=.1;
    for i to imax do
    try
     dsol(tf);
     return dsol;
    catch "cannot evaluate the solution":
            pts := Statistics:-Sample(Uniform(-eps,eps),1) ;
            ini := [ rhs(dsol(0)[2])+pts[1]] ;
            dsol('initial'=ini);
            print(i,ini)      
    end try
    end do;
    if i=imax+1 then FAIL end if;     
    end proc :

sol :=
  dsolve({ diff(u(t),t)=u(t)^2, u(0)=1},numeric,'output' = listprocedure);

sol(1);
ReInit(sol,1.0001);

sol(0);
sol(1);

You could do

(nextprime@@8)(n);

Example:

(nextprime@@8)(10000);
                             10079

Is this what you wanted?

restart;
y := 0;
for m from 0 to 2 do
   y:=`if`(m=0,x,y+(-1)^m/(m^2)^.5);
   print(m, y)
end do:

An alternative way of writing this is

restart;
y := 0;
for m from 0 to 2 do
   if m = 0 then y := x else y := y+(-1)^m/(m^2)^.5 end if;
   print(m, y)
end do:

restart;
with(LinearAlgebra):
interface(rtablesize=infinity):
node:=9:
Lambda:=DiagonalMatrix([3,6,4,2,5,2,3,5,3,7,4,7,4,6,8,9,2,1]);
Z:=Matrix(node*2):
dizLambda:=Matrix(node*2):
for i from 1 to (node*2-1) do
for j from (i+1) to (node*2) do
dizLambda[i,j]:=(Lambda[i,i]-Lambda[j,j])/2;
if dizLambda[i,j]=0 then Z[i,j]:=a elif dizLambda[i,j]<1.5 then Z[i,j]:=9 else Z[i,j]:=dizLambda[i,j] end if;
end do:
end do:

You must replace the square brackets [] used in 'de' by rounded ones ().

And since you are using numerical solution of the problem you need to give Maple the values of g, R, and omega:

g:=1: omega:=sqrt(2): R:=1:

You can use modp on the exponents of epsilon.

As an example take F to be

F:=randpoly(epsilon,degree=1000,dense):
#Then do
evalindets(F,identical(epsilon)^posint,x->applyop(rcurry(modp,2),2,x));

#You can make a procedure, which does this:
SE:=f->evalindets(f,identical(epsilon)^posint,x->applyop(rcurry(modp,2),2,x));
SE(F);

Notice that this is very fast.

Added later:

However, it seems that the following version (SE2) using 'indets' instead of 'evalindets' is faster:

SE2:=proc(f) local S1,S2,S,g;
   g:=x->applyop(rcurry(modp,2),2,x);
   S1:=convert(indets(f,identical(epsilon)^posint),list);
   S2:=map(g,S1);
   S:=zip(`=`,S1,S2);
   eval(f,S)
end proc;

You may try

F:=randpoly(epsilon,degree=10000,dense):
time(SE(F));
time(SE2(F));

Edited (once again): Reading the code (and then the help page) for 'evalindets' I found that 'evalindets' can be used with the index 'flat', which uses a faster algorithm very similar to the one I used in SE2:

evalindets[flat](F,identical(epsilon)^posint,x->applyop(rcurry(modp,2),2,x));