How to make `D` apply chain rule on user-defined symbols as it does on `Dot` and `NonCommutativeMultiply`

Maybe you can use the following:

ClearAll[dmydot, mydot]

Derivative[ndiff__][mydot][x__] := dmydot[ndiff][x];

dmydot /: Times[f__, dmydot[n1___, 1, n2___][x__]] := 
 mydot @@ MapAt[Times[f] &, {x}, Length[{n1}] + 1]

All I'm doing here is pull the inner derivative that came from the chain rule back into the mydot to replace the original entry at that position. This is done by defining the derivative as dmydot and giving it a rule that acts when something multiplies it. Now it can of course happen that this "something" also contains factors that were there even before the chain rule kicked in. But at least the operation of pulling such additional factors into the mydot product can never be incorrect. Whether you are willing to live with this possibility, I'm not sure...

Here are some tests:

D[y^2 + 2 mydot[x + b y, (x + a y), z], y]

(* ==> 2 y + 2 (mydot[b, x + a y, z] + mydot[x + b y, a, z]) *)

This now behaves exactly like Dot:

D[y^2 + 2 Dot[x + b y, (x + a y), z], y]

(* ==> 2 y + 2 (b.(x + a y).z + (x + b y).a.z) *)

Here is a more general function f inside mydot:

D[y^2 + 2 mydot[x, f[x + a y], z], y]

(* ==> 2 y + mydot[x, 2 a Derivative[1][f][x + a y], z] *)

I think this is all mathematically correct, and the main issue of unwanted inner derivatives in front of the product has been fixed. In addition, it also works with higher derivatives:

D[mydot[x, (x + a y^2), z], y, y]

(* ==> mydot[x, 2 a, z] *)

By default, Mathematica always differentiates function heads that contain the independent differentiation variables, unless the function head is listed in the "DifferentiationOptions"->"ExcludedFunctions" system options. So, we can get the behavior you want by adding mydot to this list:

SetSystemOptions[
    "DifferentiationOptions" -> 
    "ExcludedFunctions"-> DeleteDuplicates[
        Append[
            OptionValue[SystemOptions[], "DifferentiationOptions"->"ExcludedFunctions"],
            mydot
        ]
    ]
];

and then adding a custom differentiation rule for mydot:

mydot /: D[mydot[a__], y_] := D[Dot[a], y] /. Dot->mydot

Then, we have:

D[mydot[x, (x+a y), z], y]

mydot[x, a, z]

as desired. Here are the rest of @Jens' tests:

D[y^2 + 2 mydot[x + b y, (x + a y), z], y]
D[y^2 + 2 mydot[x, f[x + a y], z], y]
D[mydot[x, (x + a y^2), z], y, y]

2 y + 2 (mydot[b, x + a y, z] + mydot[x + b y, a, z])

2 y + 2 mydot[x, a f'[x + a y], z]

mydot[x, 2 a, z]