If you use square brackets around the N[i]'s then the order is respected:

solve(List,[seq(N[i],i=1..9)]);
op(op(%));
assign(%);
List;

Warnings are printed with 'print', see

showstat(WARNING);

so are not output.

Here are two ways of saving the warning messages for inspection or handling later.

The first one redefines the WARNING procedure to make it assign to the global variable 'lastwarning'.
The second method makes use of a global variable `DSI/Warnings`, which is used by 'odeplot' when _Env_in_maplet is true and is assigned the list of warnings issued. I don't use maplets so I have no idea if the assignment _Env_in_maplet:=true will have undesirable consequences in your case.

restart;
showstat(WARNING);
unprotect(WARNING);
WARNING:= proc(msg::{string, symbol}) global lastwarning;
   lastwarning:=msg;
   print(INTERFACE_WARN(1,msg,args[2 .. -1]))
end proc;
protect(WARNING);
sol:=dsolve({diff(y(x),x)=y(x)^2,y(0)=1},numeric);
plots:-odeplot(sol,[x,y(x)],-1..1);
lastwarning;
lastexception;
sol(1);
lastexception;
StringTools:-FormatMessage(lastexception[2..-1]);
WARNING(%);
lastwarning;
#Instead of redefining 'lastexception' at each call to WARNING, 'lastexception' could be assigned a cumulative list of all warnings just as is done by `dsolve/numeric/warning` in the next method: 
restart;
#showstat(`plots/odeplot`);
showstat(`dsolve/numeric/warning`);
_Env_in_maplet:=true;
sol:=dsolve({diff(y(x),x)=y(x)^2,y(0)=1},numeric);
plots:-odeplot(sol,[x,y(x)],-1..1);
`DSI/Warnings`;

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: