Thanks again, Gianluca. I am actually working on a game theoretic generalization of control theory (or Lagrangian mechanics), where instead of one functional (or running payoff) to optimize, there are two players, each with its own payoff, and each having its own "lever" or control of the common system. The Nash equilibrium at each time is complicated. There could be 0-3 solutions (the control functions A[t] and yd[t] of the two players. If it weren't for inequalities coming from physical considerations there would be 1 or 3 solutions) that are the roots of a cubic polynomial (with the coefficients of one depending on the roots of the other). An old root can go out of existence in a non-differentiable way during the time evolution of the system. Also, in some cases the boundary conditions are given at both initial and final times and new solutions can emerge even at initial times if the final time is increased.

I think the implicit differentiation method has the advantage of being "local", so that it doesn't see non-differentiable behavior until it actually happens. Perhaps it will not happen at all.

Anyway, I thought you might be interested in the background, given that you are interested in Lagrangian mechanics. I also have a "conspiracy theory" that Von Neuman was working on a formulation of Quantum mechanics that involved game theory as an alternative to the other attempts to formulate it so it gives a reasonable classical limit, given that he contributed to the foundation of both QM and game theory.

OK, maybe I will try your method since I can't figure out what is going wrong with mine. I decided to do the chain rule by hand in the differential equations. Ignore everything in comment parentheses (there are more equations that I haven't fixed yet with the new method). What I want is a way to assign Asharp=the original definition, which is a function of r[t],y[t],n[t] (as well as the other 6 dynamic variables q1r[t],q2r[t] etc) once inside the NDSolve block so it doesn't need to be computed with every reference, and also so it counts as a constant in expressions such as D[fr[r[t], y[t], n[t], Asharp], r[t]]. Same for ydsharp, Asharpdr, ydsharpdr, Asharpdn, etc, with the latter variables being computed with the implicit differentiation you suggested.

Is there a way to use Block or Module inside NDSolve to assign Asharp to have the properties I desire above?

sol2 = NDSolve[{r'[t] == fr[r[t], y[t], n[t], Asharp], 
     y'[t] == fy[r[t], y[t], n[t], Asharp, ydsharp], 
     n'[t] == fn[r[t], y[t], n[t], Asharp, ydsharp], 
     q1r'[t] == -q1r[t] D[fr[r[t], y[t], n[t], Asharp], r[t]] - 
       D[h[Asharp, r[t], y[t], n[t], ydsharp], r[t]] - 
       q1r[t] Asharpdr[
         D[[fr[r[t], y[t], n[t], A], 
            A]] /. {A -> 
             Asharp} - (Asharpdr D[h[A, r, y, n, ydsharp], 
                A] /. {A -> Asharp} + 
               ydsharpdr D[h[Asharp, r, y, n, yd], yd] /. {yd -> 
               ydsharp}), q1y'[t] == ...
         (*-q1y[t] D[fy[r[t],
         y[t],n[t],Asharp[q1r[t],q1y[t],q1n[t],q2y[t],q2n[t],r[t],y[
         t],n[t]],ydsharp[q1r[t],q1y[t],q1n[t],q2y[t],q2n[t],r[t],y[
         t],n[t]]],y[t]]-D[h[Asharp[q1r[t],q1y[t],q1n[t],q2y[t],q2n[
         t],r[t],y[t],n[t]],r[t],y[t],n[t],ydsharp[q1r[t],q1y[t],q1n[
         t],q2y[t],q2n[t],r[t],y[t],n[t]]],y[t]],q1n'[t]==-q1n[t] D[
         fn[r[t],y[t],n[t],Asharp[q1r[t],q1y[t],q1n[t],q2y[t],q2n[t],
         r[t],y[t],n[t]],ydsharp[q1r[t],q1y[t],q1n[t],q2y[t],q2n[t],r[
         t],y[t],n[t]]],n[t]]-D[h[Asharp[q1r[t],q1y[t],q1n[t],q2y[t],
         q2n[t],r[t],y[t],n[t]],r[t],y[t],n[t],ydsharp[q1r[t],q1y[t],
         q1n[t],q2y[t],q2n[t],r[t],y[t],n[t]]],n[t]],q2r'[t]==-q2r[
         t] D[fr[r[t],y[t],n[t],Asharp[q1r[t],q1y[t],q1n[t],q2y[t],
         q2n[t],r[t],y[t],n[t]]],r[t]]+D[fy[r[t],y[t],n[t],Asharp[q1r[
         t],q1y[t],q1n[t],q2y[t],q2n[t],r[t],y[t],n[t]],ydsharp[q1r[
         t],q1y[t],q1n[t],q2y[t],q2n[t],r[t],y[t],n[t]]],r[t]],q2y'[
         t]==-q2y[t] D[fy[r[t],y[t],n[t],Asharp[q1r[t],q1y[t],q1n[t],
         q2y[t],q2n[t],r[t],y[t],n[t]],ydsharp[q1r[t],q1y[t],q1n[t],
         q2y[t],q2n[t],r[t],y[t],n[t]]],y[t]]+D[fy[r[t],y[t],n[t],
         Asharp[q1r[t],q1y[t],q1n[t],q2y[t],q2n[t],r[t],y[t],n[t]],
         ydsharp[q1r[t],q1y[t],q1n[t],q2y[t],q2n[t],r[t],y[t],n[t]]],
         y[t]],q2n'[t]==-q2n[t] D[fn[r[t],y[t],n[t],Asharp[q1r[t],q1y[
         t],q1n[t],q2y[t],q2n[t],r[t],y[t],n[t]],ydsharp[q1r[t],q1y[
         t],q1n[t],q2y[t],q2n[t],r[t],y[t],n[t]]],n[t]]+D[fy[r[t],y[
         t],n[t],Asharp[q1r[t],q1y[t],q1n[t],q2y[t],q2n[t],r[t],y[t],
         n[t]],ydsharp[q1r[t],q1y[t],q1n[t],q2y[t],q2n[t],r[t],y[t],n[
         t]]],n[t]]*), r[0] == rinit, y[0] == yinit, n[0] == ninit, 
         q1r[0] == q1rinit, q1y[0] == q1yinit, q1n[0] == q1ninit, 
         q2r[0] == q2rinit, q2y[0] == q2yinit, q2n[0] == q2ninit}, {r,
   y, n, q1r, q1y, q1n, q2r, q2y, q2n}, {t, 0, 
  tfinal}(*, Method->{"Shooting", \
"StartingInitialConditions"->{r[0]==rinit,y[0]==yinit,n[0]==ninit,q1r[\
0]==q1rinit,q1y[0]==q1yinit,q1n[0]==q1ninit,q2r[0]==q2rinit,q2y[0]==\
q2yinit,q2n[0]==q2ninit}}*)]

$Aborted

It's not that I find it confusing, it's that I don't get any output. The replacements never happen, the derivative of the equations never happens, and the solutions to the equations for the derivatives in terms of the parameters never happens.In other words, the 3 commands:

equationsWithAYDdependentOnN = 
  basicEquations /. {eqs__Equal} :> ({eqs} /. {A -> A[n], 
       yd -> yd[n]});
derivativeOfTheEquationsWRTn = 
  equationsWithAYDdependentOnN /. {eqs__Equal} :> Simplify@D[{eqs}, n];
partialDerivativesOfAYD = 
 derivativeOfTheEquationsWRTn /. {eqs__Equal} :> 
   FullSimplify@First@Solve[{eqs}, {A'[n], yd'[n]}]

do nothing.

Yes, you want the delayed substitution. With immediate substitution the operations are made before the pattern-matching is made, which in your case does not make sense, because you replace the variables A, yd inside the symbol eqs, which achieves nothing at all. You must first find the equations matching the pattern, and then make the replacements.

If you find all this too confusing, you may look into the Piecewise expression and copy-and-paste the equations manually into other cells, where you can make the replacements. You can keep track of the Piecewise conditions manually too. I suggest to make the algebraic calculations with exact rational coefficients, because Solve feels more comfortable that way.

With one underscore it matches only a single equation.

In your case there are more than one equations, and you need two underscores.

Yes, the eqs__Equal is a pattern that isolates the equations inside the Piecewise structure. When you take the derivative of an equation, it automatically takes the derivatives of the two sides. This derivative is itself an equation that contains A'[n] linearly. You solve the equation for A'[n]. You get a formula that contains A[n].

Actually, your original equations in A,nd are not hopelessly hard to solve symbolically. You may give it a try:

basicEquations /. {eqs__Equal} :> Solve[{eqs}, {A, yd}]

If all your variables are symbolic you get very complicated algebraic formulas, but for specific numeric values of the parameters they may become manegeable.

I tried copying your code, but I don't get derivatives at all. I get equations for A and yd (in terms of each other, so still implicit). What is "eqs_Equal"? I have an older version and maybe this is a new thing.

Thanks, but I don't know if this applies to my situation. I don't think it is possible to solve for the derivatives analytically, as in this example, but I could be wrong. Here is the code:

fr[r, y, n, A] := A (rcg[hi[r, y, n, A]] - krg3 ncd[rcg[hi[r, y, n, A]]]) + (1 - A) rig[y] - A rcd - (1 - A) rid fy[r, y, n, A, ycd_] := r p[r, y, n, A] (A ycg[rcg[hi[r, y, n, A]], ycd] + (1 - A) yig[n]) fn[r, y, n, A, ycd_] := nig + ncg* n (h[A, r, y, n, ycd] - hc[ycd]) - (A ncd[rcg[hi[r, y, n, A]]] + nid)

rcg[h_] := krg1 - (h krg2 )
ycg[rp_, yd_] := kyg1 - rp kyg2 - kyg3 yd
ncd[rp_] := knd1 rp^2 
hi[r_, y_, n_, 
  A_] := (kr1 r - kr2 r^2) (ky1 y - ky2 y^2 ) (kn1 n - kn2 n^2) (1 + 
    khi A)
h[A_, r_, y_, n_, yd_] := A hc[yd] + (1 - A) hi[r, y, n, A]
rig[y_] := kig1 + y kig2 
hc[yd_] := hci - kc yd^2 
p[r_, y_, n_, 
  A_] := (hi[ky1/(2 ky2), kr1/(2 kr2), kn1/(2 kn2), A] - 
    hi[r, y, n, A] )/(hi[ky1/(2 ky2), kr1/(2 kr2), kn1/(2 kn2), A] - 
    hi[rinit, yinit, ninit, A])
yig[n_] := kyig1 + n kyig2 
ydsharpquad[q2y_, q2n_, r_, y_, n_, A_] :=

 Min[-0.01, 
  Coefficient[(q2y - 1) fy[r, y, n, A, w] + q2n fn[r, y, n, A, w](*,
   hc[w]>=0*), w^2]]

ydsharplin[q2y_, q2n_, r_, y_, n_, A_] :=

 Coefficient[(q2y - 1) fy[r, y, n, A, w] + q2n fn[r, y, n, A, w](*,hc[
  w]>=0*), w] 
ydtent[q2y_, q2n_, r_, y_, n_, 
  A_] := -0.5 ydsharplin[q2y, q2n, r, y, n, A]/
   ydsharpquad[q2y, q2n, r, y, n, A]
Asharpquad[q1r_, q1y_, q1n_, r_, y_, n_, yd_] := 
 Min[-0.01, 
  Coefficient[
   q1r fr[r, y, n, w] + q1y fy[r, y, n, w, yd] + 
    q1n fn[r, y, n, w, yd] + h[w, r, y, n, yd](*,0<=w<=1*), w^2]]
Asharplin[q1r_, q1y_, q1n_, r_, y_, n_, yd_] := 
 Coefficient[
  q1r fr[r, y, n, w] + q1y fy[r, y, n, w, yd] + 
   q1n fn[r, y, n, w, yd] + h[w, r, y, n, yd](*,0<=w<=1*), w]
Atent[q1r_, q1y_, q1n_, r_, y_, n_, 
  yd_] := -0.5 Asharplin[q1r, q1y, q1n, r, y, n, yd]/
   Asharpquad[q1r, q1y, q1n, r, y, n, yd]
sol[q1r_, q1y_, q1n_, q2y_, q2n_, r_, y_, n_] := 
 Solve[{A == Atent[q1r, q1y, q1n, r, y, n, yd], 
   yd == ydtent[q2y, q2n, r, y, n, A], 0 <= A <= 1, yd >= 0}, {A, yd},
   Reals]
Asharp[q1r_, q1y_, q1n_, q2y_, q2n_, r_, y_, n_] := 
 If[sol[q1r, q1y, q1n, q2y, q2n, r, y, n] == {}, 0, 
  sol[q1r, q1y, q1n, q2y, q2n, r, y, n][[1]][[1]][[2]]]
ydsharp[q1r_, q1y_, q1n_, q2y_, q2n_, r_, y_, n_] := 
 If[sol[q1r, q1y, q1n, q2y, q2n, r, y, n] == {}, 0, 
  sol[q1r, q1y, q1n, q2y, q2n, r, y, n][[1]][[2]][[2]]]
q1rinit = 0.1
q1yinit = 0.1
q1ninit = 0.1
q2rinit = -0.1
q2yinit = -0.1
q2ninit = -0.1
krg1 = 1; krg2 = 0.1; kyg1 = 1; kyg2 = 0.1; kyg3 = 0.1; knd1 = 1; kr1 \
= 1; kr2 = 0.01; ky1 = 1; ky2 = 0.01; kn1 = 1; kn2 = 0.01; kig1 = 1; \
kig2 = 1; hci = 1.1; krg3 = 1; rcd = 1; rid = 0.1; ncg = 1; nig = 1; \
kc = 0.1; yinit = 1; rinit = 1; ninit = 1.; nid = 1; tfinal = 0.01; \
kyig1 = 1; kyig2 = 0.1; yid = 1; khi = 0.5; Ainit = 0.5; ydinit = 0; \
Afinal = 0; ydfinal = 0.1

0.1

0.1

0.1

-0.1

-0.1

-0.1

0.1
In[391]:= 
Derivative[0, 0, 0, 0, 0, 0, 0, 1][
  Asharp][q1rinit, q1yinit, q1ninit, q2yinit, q2ninit, rinit, yinit, \
ninit]

Out[391]= $Aborted
In[398]:= Asharp[q1rinit, q1yinit, q1ninit, q2yinit, q2ninit, rinit, \
yinit, ninit]

During evaluation of In[398]:= Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

Out[398]= 0.175365