So I am solving a PDE for a function $h(\theta,t)$ via finite difference scheme. The PDE has a function $Q$ in it, which I wish to update each time step depending on where $h$ lies.

Firstly, let us look at simple case when $Q = const$. Here is my code:

```
ClearAll[dh, t, n, ds, h0, sol, h, valh, ti, s];

n = 1001;
ds = Pi/(n - 1);

(* dh: R^n \[Rule] R^n, symbolic and numeric *)
dh[h_List /; Length@h == n] := With[{

    s = Array[# &, n, {-π, 0}],
    Q = 0.005,(*I wish to update this each time step depending on \
where h lies.*)
    d1 = ListCorrelate[{-0.5, 0, 0.5}/ds, #] &,
    d2 = ListCorrelate[{1, -2, 1}/ds^2, #] &,
    h2 = Join[{2}, h, {2}]
    }, (-2 + h) (-1 + 
      h) (-18 Q d1[h2] + (-2 + h) (-1 + h) (4 + h (-15 + 8 h)) Cos[
        s] d1[h2] + (-2 + h) (-1 + 
         h) ((-9 + 6 h) d1[h2]^2 - (2 - 3 h + h^2) (h - d2[h2])) Sin[
        s])
   ];

(*h0 ∈ R^n *)
h0 = Table[Piecewise[{

     {1.999, -Pi < θ < -7 Pi/8},
     {(2 Sin[-6 Pi/7])/
       Sin[θ], -6 Pi/
         7 < θ < -2.0909003694838066`}, {1.001, \
-2.26998149133592` < θ < -0.6662394324925154`},
     {(2 Sin[-(Pi/10)])/
       Sin[θ], -0.6662394324925154` < θ < -(Pi/10)},
     {1.999, -(Pi/8) < θ < 0}

     }, 1.999]

   , {θ, -1. Pi, 0., ds}
   ];

(*Solve 
dh/dt(t) = dh(t)
h(0)= h0

for h(t) ∈ R^n
for t ∈ (0,20)
*)
sol = With[{
    h = h[#][t] & /@ Range@n
    },

   NDSolveValue[
    {
     D[h, t] == dh[h],
     (h /. t -> 0) == h0
     },
    Head /@ h,
    {t, 0, 20}
    ]

   ];

valh = Transpose[#["ValuesOnGrid"] & /@ sol];
ti = Flatten@sol[[1]]["Grid"];
s = Array[# &, n, {-Pi, 0}];

h = ListInterpolation[valh, {ti, s}, InterpolationOrder -> 1];

ListAnimate@Table[
  Show[
   (*Interface*)
   PolarPlot[
    Evaluate[h[ti, θ]],
    {θ, -π, 0},

    PlotRange -> {{-2, 2}, {-2, 2}},
    ImageSize -> Large,
    PlotStyle -> Black],
   (*Ring*)
   RegionPlot[
    1 <= x^2 + y^2 <= 4,
    {x, -2, 2}, {y, -2, 2},

    PlotStyle -> Opacity[0.1], AspectRatio -> 1
    ]

   ]
  , {ti, 0, 20, 0.2}]

```

[![enter image description here](https://i.stack.imgur.com/gdjtR.gif)](https://i.stack.imgur.com/gdjtR.gif)

As we can see the interface moves around with a fixed motion, due to Q being fixed at 0.005. However, in reality I wish to update $Q$ at each time step with the following relation $Q = \cos \psi - \cos \phi$, where $\phi$ and $\psi$ are the angles from the base of the annulus to the middle of the right and left flat sections.

[![enter image description here](https://i.stack.imgur.com/T7Mhw.png)](https://i.stack.imgur.com/T7Mhw.png)

Now I use the following pattern matching code to find the initial $Q$. It basically looks for the median position of the centre of the interface and multiply by `ds` to convert it to an angle.

```
Q = N@Cos[Median@Position[h0[[1 ;; 500]], x_ /; 1.001 < x < 1.999] ds][[1]] - N@Cos[(Median@
  Position[h0[[501 ;; 1001]], x_ /; 1.001 < x < 1.999]) ds][[1]] 

```

This gives me a sensible initial $Q$ of 0.26. I wish to then update this $Q$ with each time step, with the new position of the middle of the interface. My hope is eventually that the two angles will become equal and $Q = 0$. Here is my attempt to update $Q$ at each time step in the Finite Difference scheme.

```
n = 1001;
ds = Pi/(n - 1);
ClearAll[dh, t];
dh[h_List] := 
  With[{s = Array[# &, n, {-π, 0}],
    Q = N@Cos[Median@Position[#[[1 ;; 500]], x_ /; 1.001 < x < 1.999] ds][[1]]-N@Cos[(Median@Position[#[[501 ;; 1001]], x_ /; 1.001 < x < 1.999]) ds[[1]] &,
    d1 = ListCorrelate[{-0.5, 0, 0.5}/ds, #] &,  
    d2 = ListCorrelate[{1, -2, 1}/ds^2, #] &},
    (-2 + h) (-1 + h) (-18 Q[#1] d1[#1] + (-2 + h) (-1 + h) (4 + h (-15 + 8 h)) Cos[
      s] d1[#1] + (-2 + h) (-1 + 
       h) ((-9 + 6 h) d1[#1]^2 - (2 - 3 h + h^2) (h - 
          d2[#1])) Sin[s]) &@Join[{2}, h, {2}]];

```

Which runs, but nothing moves! It should really as the initial $Q$ for $h\_0$ is 0.26 as shown before. 

Any pointers on how to update $Q$ at each time step? I feel like I'm nearly there but perhaps my syntax is incorrect? Maybe using the pure function incorrectly?

Any help appreciated as always.

**Edit**

So another idea I had was to include $Q$ as an unknown in the finite difference scheme. ie. It would appear as `Q[1]`, `Q[2]` etc. in each of the equations. I could then solve the finite difference scheme `D[h,t] == dh[h,Q]` coupled with an equation for $Q$, in the `NDSolveValue`. Here is my attempt:

```
n = 1001;
ds = Pi/(n - 1);
ClearAll[dh, t];
dh[h_List, Q_] := 
  With[{s = Array[# &, n, {-π, 0}],
     d1 = ListCorrelate[{-0.5, 0, 0.5}/ds, #] &, 
     d2 = ListCorrelate[{1, -2, 1}/ds^2, #] &},
(-2 + h) (-1 + 
     h) (-18 Q d1[#1] + (-2 + h) (-1 + h) (4 + h (-15 + 8 h)) Cos[
       s] d1[#1] + (-2 + h) (-1 + 
        h) ((-9 + 6 h) d1[#1]^2 - (2 - 3 h + h^2) (h - 
           d2[#1])) Sin[s]) &@Join[{2}, h, {2}]];

q0 = 0.26

Clear@h;
sol = With[{h = h[#][t] & /@ Range@n, Q = Q[#][t] & /@ Range@n}, 
NDSolveValue[{D[h, t] == dh[h, Q], 
  Q[t] == -N@
      Cos[(Median@
           Position[h[[501 ;; 1001]], 
            x_ /; 1.001 < x < 1.999]) ds][[1]] + 
    N@Cos[Median@
         Position[h[[1 ;; 500]], 
          x_ /; 1.001 < x < 1.999] ds][[1]] , (h /. t -> 0) == 
   h0, (Q /. t -> 0) == q0}, Head /@ h, {t, 0, 20}, 
 Method -> {"EquationSimplification" -> "Residual"}]];

```

But this gives me a very strange error.

> 
> `Transpose::nmtx`: The first two levels of `{NDSolve`xs$422676, NDSolve`xs$422677, NDSolve`xs$422678, NDSolve`xs$422679, NDSolve`xs$422680, NDSolve`xs$422681,NDSolve`xs$422682, NDSolve`xs$422683, NDSolve`xs$422684, etc. <<951>>}` cannot be transposed.
> 
> 
>