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;

You could use translate from the plottools package or the more general transform from the same package:

restart;
with(plottools);
with(plots):
d1:=disk([1,2],1,color=blue):
display(d1);
a:=2.52:
display(d1,translate(d1,-a,3*a),scaling=constrained);
#Alternatively use the more general transform:
tr:=transform((x,y)->[x-a,y+3*a]):
display(d1,tr(d1),scaling=constrained);

indets(sys);
A,B:=LinearAlgebra:-GenerateMatrix(sys,[a, b, c, d, e, g]);