This skeleton example shows how to do what you want.  You can jazz it up as needed.

Download mw.mw

I have made many corrections and changes to your worksheet.  This may or may not be what you really need, therefore look at everything closely!

Download mw.mw

Your f is a function of two variables.  What do you mean by "the derivative of f"?

The partial derivative of f with respect to its first variable (which is t) is written D[1](f).  To evaluate that derivative at the point (t[n-1],w[n-1]), we write
D[1](f)(t[n-1],w[n-1]).

If you need the derivative of f with respect to its second variable (which is y), then you will need
D[2](f)(t[n-1],w[n-1]).

I don't know which of two you need in your code because it's not clear to me what it is attempting to do.

In the first line of the image that you have attached, we see the input is
f := whatever
but the corresponding output is
f = whatever

Note the disagreement between the ":=" and "=" in the input versus the output.  This means that you have first executed the input without a colon and obtained the output, and then you have inserted the colon in the input but have not executed it!  As a result, f has not been assigned a value and therefore it appears as a plain f in the final result.

Execute the worksheet again without changing anything.  You should get the expected result.

Instead of plotting and then shifting, you may want to plot the arrows directly in the desired location.

For instance, to plot arrows originating at the point [3,4], you can do:

restart;
plots:-arrow([3,4], [[1,0], [0,1]]);

Moreover, as we see here, a single call to arrow is sufficient to plot any number of arrows.

With the unknown function u(x,y), the Dirichlet condition  u(0,y) = f(y) is entered as

u(0,y) = f(y);

The Neumann condition ∂u/∂x(0,y) = g(y) is entered as 

D[1](u)(0,y) = g(y);

The Neumann condition ∂u/∂y(x,0) = h(x) is entered as 

D[2](u)(x,0) = h(x);

The notations D[1] and D[2] refer to differentiation with respect to the first and second variables of u(x,y).

The presence of the cosine term is more of an effect of the geometry than physics, although I dislike putting a hard barrier between the two subjects.

The Pythagorean Theorem says that the square of the hypotenuse is the sum of the squares of the two other sides in a right triangle.  The generalization to an arbitrary triangle is known as the law of cosines.  I trust that you already know that but in the unlikely case, look it up in Wikipedia.

The law of cosines enters your calculations because the two vectors that you have shown are not perpendicular. The pearl's velocity is the sum of those two vectors.  Calculating the sum relies on the law of cosines.

Download mw.mw

This works in Maple 2019.  It may require minor adjustments in older versions.

Download mw.mw

As tomleslie said, there are many ways of achieving the result.  Here is a straightforwad "bare hands" way which may be more accessible to a beginner.

Download mw.mw

Probably this is what you want:

for n in (10*2^i $i=0..4) do
    print(n);  # or whatever
end do;

or equivalently

for i from 0 to 4 do
    n := 10*2^i;
    print(n); # or whatever
end do;

or

n := 10; 
while n <= 160 do
    print(n); # or whatever
    n := 2*n;
end do;

Pick whichever you like.

Here is a possible way.

Download lhs_rhs-new.mw

I have made several changes to your worksheet and I hope that I have not introduced any errors.  Check to be sure.

The optimal diffusion coefficient turns out to be d = 0.14 (approximately).

Download zz.mw

To plot the derivatives, add the output=operator option to your dsolve(), as in
p := dsolve(%, vars, numeric, output = operator);

Do odeplot() as before to plot the solutions.  To plot the derivative of y(t) (also those of x[1](t) and x[2](t)) do
eval(diff(y(t), t), Sys);
eval(%, p);
plot(%, t = 0 .. 1500);

As to the higher order derivatives, y'', y''', etc., these are not quite well-defined because the right-hand sides of your differential equations contain discontinuities, which means that y' is discontinuous.  You don't want to differentiate a discontinuous function, do you?

You can actually see the discontinuity in y' by changing the plot(%, t = 0 .. 1500) in what I have shown to plot(%, t = 500 .. 1500).

As to the plotting of the bar graph, I don't think it's difficult but I don't have the time to go into it right now.  Perhaps someone else will.

Of the two solutions returned by dsolve() one corresponds to positive y and the other to negative y.

restart;
de := diff(y(x),x) = abs(y(x)) + 1;
                         d                   
                  de := --- y(x) = |y(x)| + 1
                         dx                  
dsol := dsolve(de);
                       exp(-x)                           
      dsol := y(x) = - ------- + 1, y(x) = exp(x) _C1 - 1
                         _C1                             
odetest(dsol[1], de) assuming y(x)<0;
                               0
odetest(dsol[2], de) assuming y(x)>0;
                               0

Download mw.mw