Array operations, what's allowed and what's not?

[JianghuiDu]

I can't find any documentation on non-scalar usage. Only after reading some issues I got to know a bit more about array operations. But still there are many unclear things.
using a simple example:

Ngrid = 100
# A_POC, A_O2 are constant tridiagonal matrices...
@variables t POC[1:Ngrid](t) O2[1:Ngrid](t) RO2POC[Ngrid](t)
@parameters p1, p2
D = Differential(t)
f! = function(x) x end

What I want is something like

@named modeltest = ODESystem(
    [
    mul!(D.(POC),A_POC,POC);
    mul!(D.(O2),A_O2,O2);
    D(POC[1]) ~ D(POC[1]) + POC[1]/10 + 1.0;
    RO2POC .~ O2 ./ (100.+ O2) .* p1 ./ (p2 .+ 10) .* POC;
    D.(POC) .~  D.(POC) .- RO2POC;
    D.(O2) .~ D.(O2) .- RO2POC;
    f!.(D.(O2))
]
)

1. It looks like, mul!(D.(POC),A_POC,POC) is not allowed? D.(POC) .~ mul!(A_POC,POC) is also not allowed? Seems one must do D.(POC) .~ A_POC*POC? (also I didn't see any mention of .~ or D.() in the docs.). Does not using mul! degrade performance?
2. The derivatives can only be assigned once and no modification allowed later? Like after D.(POC) .~ A_POC*POC one cannot do D.(POC) .~ D.(POC) .-RO2POC?
3. It is said that the temporary variable like RO2POC will be tracked in the modelling process even that it will be eliminated in the final model code. Does that consume more memory and reduce performance. Can that be disabled?

[ChrisRackauckas]

These are being changed by array symbolics so it'll be fairly different.