You have asked how to calculate the partial derivative of u with respect to y.  Here is how.

Download numerical-flux.mw

Disclaimer: This is not an answer!

I am posting this as an "Answer" to set it apart from the twisty discussion thread that we have had so far.

I have converted the OP's equations and data into a self-contained Maple worksheet.  At the end of the worksheet I have inserted the questions that the OP is seeking answers for.  My attempts to solve the equations were unsuccessful.  See what you can do with this.

Download rotor-equation.mw

What you have shown appears to be an incomplete fragment (what is "Thust", for instance?) and therefore I cannot tell you how to fix that.  If Maple's fsolve cannot solve your equation, I doubt that your own Newton's root-finder will do any better.

Nonetheless, since you asked for a recursive implementation of Newton's root-finding method, here it is, for whatever it's worth.

Download newton-iteration.mw

Download mw.mw

TomLeslie has already answered your question. Here is a shorter alternative.

Download mw.mw

I don't know why you refer to singular solutions. I don't see anything "singular" there.

Neither of the two solutions produced by your construction is a solutions.  You are "misleading" Maple into looking at the ODE as one of the d'Alembert type.  Maple blindly obliges and produces the wrong results.

Here is how one (illegally) casts your equation into a d'Alembert type.

1. Square both sides of the equation to get (dy/dx)^2 = 1 + x + y;
2. Rearrange into y + 1 = - x + (dy/dx)^2
3. Let z = 1 + y.  Then z = - x + (dz/dx)^2

This is a d'Alembert ODE (see ?Solving d'Alembert ODEs) with f(s)=-1 and g(s)=s^2.

Where is the error?  It is in Step 1.  By squaring the two sides of the ODE, we are introducing spurious solutions.  Maple's solution would be correct (for the wrong reasons) if you replace the right-hand side of your ODE by its negative.

To avoid the issue, don't impose the dalembert requirement.  Do dsolve(ode_1) to let Maple do its own thing and produce the correct, albeit implicit, solution.

This comes close to what you have sketched.

p := plots:-display(
	plot(floor(x), x=1..8, discont=[color="Red"], color="Green"),
	plot(floor(x), x=-8..-1/2, discont=[color="Red"], color="Green"),
	symbol=solidcircle, symbolsize=10, thickness=3):
plottools:-transform((x,y)->[1/x, 1/y])(p);

Download mw.mw

Your f,fdiff=... command is equivalent to the simpler:

f, diffs := GenerateMatrix(sys222, var1);

However...

A matrix representation of a system of equations is useful if the system is linear. Your system is nonlinear. What do you expect to get out of a matrix representation of a nonlinear system?

Regarding the symbolic option to simplify(), the help page says:

The result of such an operation is in general not valid over the whole complex plane and can lead to incorrect results if you assume the expressions represent analytical functions.

Thus, the symbolic option should be used only when you really know what you are doing.  Preferably it should not be used at all.  In your worksheet you are comparing two simplifications, one with and one without the symbolic option.  It's not surprising that the result are inconsistent.

restart;
doit := proc(oper)
  local x := 3, y := 4;
  return oper(x,y);
end proc:
doit(`+`);
                               7
doit(`-`);
                               -1
doit(`*`);
                               12
doit(`/`);
                               3/4

The error message says it all; it doesn't know what to do with d.

In the IBC you have

C(0, t) = 1 + d(1 - cos(varpi*t))

What is d?  If it's a function, then you need to specify it.  If it's a multiplier constant, then you need to provide its value and also don't forget to insert a multiplication sign between d and the adjacent parenthesis.

Regarding the differential equation for theta:

It has A[1] and A[2] in it, but there are no A[1] and A[2] in your expected solution.  Is that an oversight?

By the way, that differential equation is quite trivial to solve by hand.  It's not worth bothering with Maple with it. If you try the manual solution, you will find that a solution exists provided that  - v[1]  - A[1] + A[2] = v[3]. 

Regarding the differential equation for phi:

Express the differential equation in terms of Heaviside() rather than pieceswise(), as in:
 

eq2 := diff(phi(t),t) = -v[1] + v[1]*Heaviside(t-T[1]) + v[3]*Heaviside(t-T[2]);
dsolve({eq2,phi(0)=-v[2]});

That yields the unique solution of that initial value problem, and therefore determines the value of phi(1).  You cannot arbitrarily demand a different value for phi(1).
 

As dharr has correctly pointed out, a nonzero initial velocity obviates the objection that I had expressed in my prevrious message.  Following through with that insight, we end up with a semi-cubic parabolic curve y^3 = - a*x^2 for the particle's path.  The adjustment, relative to dharr's solution, is that the x compoment of the initial vecocity should be zero rather than v₀.

Download semi-cubic-path.mw