Try this instead.

restart;
with(plots):
MW := 29; #assignment uses :=
cp := 1200;
hf := 4*10^7;
AF := 16;
eq1 := (D(F))(t) = 6.19*10^9*exp(-15098/T(t))*F(t)^.1*(.2*O(t))^1.65; #T(t)!!!! not T
eq2 := (D(O))(t) = 16*(D(F))(t);
eq3 := (D(Pr))(t) = -17*(D(F))(t);
eq4 := (D(T))(t) = (D(F))(t)*hf*(8314*T(t))/((cp-8314)*P(t));
eq5 := (D(P))(t) = P(t)*(D(T))(t)/T(t)-P(0)*(D(T))(t)/T(0);
ics:={O(0) = 152.32*10^(-3), P(0) = 101325, T(0) = 753, F(0) = 9.52*10^(-3), Pr(0) = 0};
#Equations eq2 and eq3 ought to be left out since O(t) and Pr(t) are trivially found from F(t):
#eq2:
eval(O(0)=16*F(0)+C,ics);
solve(%,C);
#eq3
eval(Pr(0)=-17*F(0)+C,ics);
c:=solve(%,C);
#Thus
OF:=O(t)=16*F(t);
PrF:=Pr(t)=-17*F(t)+c;
#The reduced system:
ics2:={P(0) = 101325, T(0) = 753, F(0) = 9.52*10^(-3)};
sys2:=eval({eq1, eq4, eq5},{OF,PrF} union ics2);
sol := dsolve(sys2 union ics2, numeric);
odeplot(sol,[[t,P(t)]],0..10^6);
odeplot(sol,[[t,F(t)]],0..10^6);
odeplot(sol,[[t,T(t)]],0..10^6);

A:=Matrix([[8,4,2,1],[4,8,2,1],[2,2,8,1],[1,1,1,8]]);
with(LinearAlgebra):
IsDefinite(A);
L,U:=LUDecomposition(A,output=['L','U']);
D1:=DiagonalMatrix([seq(U[i,i],i=1..4)]);
L.D1.Transpose(L);

After differentiating try
simplify(%);
convert(%,cot);

You could use LeastSquares from the LinearAlgebra package.
If the entries are rational numbers then use the optional argument 'optimize' (or optimize=true). Are they floats this seems unnecessary.
optimize=true minimizes the 2-norm among the infinitely many solutions (see the help page for LinearAlgebra:-LeastSquares).

restart;
with(LinearAlgebra):
A:=RandomMatrix(5,10,outputoptions=[datatype=float]);
b:=RandomVector(5,outputoptions=[datatype=float]);
LinearAlgebra:-LinearSolve(A,b);
LinearAlgebra:-LeastSquares(A,b);
X:=LinearAlgebra:-LeastSquares(A,b,optimize);
Arat:=Matrix(convert~(A,rational));
brat:=Vector(convert~(b,rational));
LinearAlgebra:-LeastSquares(Arat,brat);
Xr:=LinearAlgebra:-LeastSquares(Arat,brat,optimize);
evalf[15](Xr)-X;

Define A to be the matrix
A:=Matrix([[0,1],[-3,c]]);
then your system can be written x' = A.x, where x is the vector <x1,x2>.
This linear system has <0,0> as its equilibrium point. So the question is about the change in qualitative behavior as c changes: Stability and type of point.
To answer this you have to look at
det:=Determinant(A);
Tr:=Trace(A);
dscr:=Tr^2-4*det;
And surely this must ring a bell.

restart;
#Introduce b=c/m for simplicity:
sys:={diff(x(t),t,t)=-b*sqrt(diff(x(t),t)^2+diff(y(t),t)^2)*diff(x(t),t),
   diff(y(t),t,t)=-b*sqrt(diff(x(t),t)^2+diff(y(t),t)^2)*diff(y(t),t)};
ics:={x(0)=0, y(0)=0, D(y)(0)=V0*sin(theta), D(x)(0)= V0*cos(theta)};
dsolve(sys union ics); #Direct attack doesn't work
#Introduce the speed v(t):
eq:=v(t)=sqrt(diff(x(t),t)^2+diff(y(t),t)^2);
diff(eq,t);
eval(%,sys);
simplify(%);
#So the speed satisfies this simple first order ode:
ode_v:=subs(diff(x(t),t)^2+diff(y(t),t)^2=v(t)^2,%);
dsolve({ode_v,v(0)=V0}); #The speed has been found
subs(rhs(eq)=rhs(%),sys); #Insert result in the original system
dsolve(% union ics); #Solution found
plot(eval([cos(theta)*ln(V0*b*t+1)/b,sin(theta)*ln(V0*b*t+1)/b],{V0=1,b=1,theta=Pi/3}),t=0..10);
#Can of course be done numerically too, but then you need values for the parameters.
#These can be changed, though, via the parameters option:
res:=dsolve(sys union ics,numeric, parameters=[b,V0,theta],output=listprocedure);
X,Y,Xp,Yp:=op(subs(res,[x(t),y(t),diff(x(t),t),diff(y(t),t)])); #Not needed here
res(parameters=[1,1,Pi/3]): #Setting parameters
plots:-odeplot(res,[[t,x(t)],[t,y(t)]],0..10,thickness=2);

I changed the definitions of parDist and perpDist to the following piecewise versions:

parDist:=(x,y)-> piecewise( evalf(y)<evalf(-tan(Pi/90)*(x-centralR)),evalf(centralR*(Pi/12+arctan(y/x))),evalf(sin(Pi/90)*(x-centralR)+cos(Pi/90)*y+centralR*(arctan(y/x)+Pi/12)));
perpDist:=(x,y)-> piecewise(evalf(y)<evalf(-tan(Pi/90)*(x-centralR)),evalf(sqrt(x^(2)+y^(2))-centralR),evalf(cos(Pi/90)*(x-centralR)-sin(Pi/90)*y));


and tried
FourierSeries(ByFin(x, t, 0), [t = 0.1e-1 .. 1, x = 1.5 .. 1.6], [5, 5], output = numeric, series = sine);
as well as
FourierSeries(ByFin(x, t, 0), [t = 0.1e-1 .. 1, x = 1.5 .. 1.6], [5, 5], output = numeric, series = cosine);

Both worked fast.
For your convenience I upload the altered version of your worksheet:

MaplePrimes13-1019Or.mw

K:=sin(a+b)^2-sin(a-b)^2; #Surely you didn't use the syntax you have in your question?
selectremove(has,expand(K),a);
combine~([%]);
`*`(op(%));

You may try the link
http://www.mapleprimes.com/questions/81060-Muller-Method

In view of the RootOf results of these very special cases, I rather doubt it:

solve(eval(f,{n=1}),beta[1]);
solve(eval(f,{n=2}),beta[1]);
solve(eval(f,{n=1,theta[1]=2}),beta[1]);

It is easier to use DeleteColumn (and your vector b is wrong):
b:=Column(A,3); #You had 2
C:=DeleteColumn(A,3);
X:=LinearSolve(C,b);
X[8];

Since you set z:=ic;  the loop doesn't run even once because of the while clause.
If the intention is to skip the while check at the first iteration you could do
for t from 1 to Nit while z<>ic or t=1 do
etc.
I included ... to Nit... to prevent the loop from running forever.
Alternatively, you could put an extra command
z:=Map(z);
before the loop starts. But I would still want to make sure that the loop doesn't run forever.
I didn't try to check the code, primarily because I could not copy your text since it is presented as an image.

Here is a way to get at least part of what you want and in some respects more than you want.
I start by solving the problem for a single initial speed. After that I make a procedure that takes an initial speed and outputs the position vector corresponding to that speed.
I modified your Eskritt to do the system as a whole, i.e.
dv/dt = F/m
dr/dt = v
No record is kept of the velocity vector.

restart;
C_0:=0.4: d:=0.0002: mu:=1.8*10^(-5): rho:=1.2: g:= 9.81:
rho_v:=1000.0: m:=evalf((4/3)*Pi*rho_v*(d/2)^3):

F_D:=proc(v::list) # kraft/masse
  local k,vx,vy;
  vx:=v[1];
  vy:=v[2];
  k:=-evalf((3*Pi*mu*d+(Pi/8)*C_0*rho*(d^2)*sqrt(vx^2+vy^2))/m);
  [k*vx,k*vy]
end proc;

Eskritt:=proc(vo::list,h)
  local vx,vy, F;
  F:=F_D(vo);
  vx:=vo[1]+h*F[1];
  vy:=vo[2]+h*(F[2]-g);
  [vx,vy]
end proc;
phi:=Pi/5: v0:=evalf([v_0*cos(phi), v_0*sin(phi)]): v[0]:=v0;
h:=0.002: N:=150:
##################################
#First we take a single initial speed for illustrative purposes:
r:=Vector(N): #r will be a vector having positions (lists) as entries
v_0:=4: #Initial speed (example)
v1:=v0: #Initial velocity (list)
r[1]:=[0,0]: #Initial position
#The while clause is used to determine when the ball hits the ground:
for i from 2 to N while r[i-1][2]>=0 do
  r[i]:=r[i-1]+h*v1;
  v1:=Eskritt(v1,h); #No record is kept of the velocity
end do:
i;
r[i-1];
r[i-2];
plot(r[1..-3],labels=[x,y],caption=typeset("Initial speed = ",v_0));
#################################
#Now the procedure mentioned in the beginning:
R:=proc(speed) local r,v1,i; global v_0,v0,N,h;
   r:=Vector(N); #r will be a vector having positions (lists) as entries
   v_0:=speed; #Initial speed
   v1:=v0; #Initial velocity (list)
   r[1]:=[0,0]; #Initial position
   for i from 2 to N while r[i-1][2]>=0 do
      r[i]:=r[i-1]+h*v1;
      v1:=Eskritt(v1,h)
   end do;
   if i>N then r else r[1..i-1] end if
end proc;
#Test:
R(4);
plot(R(4));
#An animation with traces:
plots:-animate(plot,['R(v)'],v=0.2..4,frames=20,trace=20);

From the error you notice that h and v are not given values.
You do assign to h inside the procedure wave, but you are  assigning to the local h, and it is the global h that appears in phi.
You could remove h from the locals list and declare it global:
global h;
However, inside wave you have the expression
phi(u_past[j-1],u_past[j],u_past[j+1],s,h)
which has h as the last argument to phi, yet you have defined phi globally without referring to h:
phi:=unapply(rhs(phi),u_S,u_W,u_E,i,j);

What is v? That appears in phi as well.

You should use definite integrals:
kinematics := proc(a::procedure, t::name, s0::numeric, v0::numeric, t0::numeric, t1:: numeric)
  local tau,v:=t->v0+int(a(tau),tau=t0..t),s:=t->s0+int(v(tau),tau=t0..t);
  plot([a(t),v(t),s(t)], t=t0..t1);
end proc;