Update docs

Author: YingboMa

docs/src/highlevel.md

@@ -10,7 +10,7 @@ methods.
 ```@docs
 @parameters
 @variables
-@derivatives
+Differential
 Base.:~(::Num, ::Num)
 modelingtoolkitize
 ```

docs/src/tutorials/higher_order.md

@@ -15,7 +15,7 @@ using ModelingToolkit, OrdinaryDiffEq
 
 @parameters t σ ρ β
 @variables x(t) y(t) z(t)
-@derivatives D'~t
+D = Differential(t)
 
 eqs = [D(D(x)) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
@@ -25,10 +25,11 @@ sys = ODESystem(eqs)
 ```
 
 Note that we could've used an alternative syntax for 2nd order, i.e.
-`@derivatives E''~t` and then `E(x)` would be the second derivative,
-and this syntax extends to Nth order.
+`D = Differential(t)^2` and then `E(x)` would be the second derivative,
+and this syntax extends to `N`-th order. Also, we can use `*` or `∘` to compose
+`Differential`s, like `Differential(t) * Differential(x)`.
 
-Now let's transform this into the ODESystem of first order components.
+Now let's transform this into the `ODESystem` of first order components.
 We do this by simply calling `ode_order_lowering`:
 
 ```julia

docs/src/tutorials/ode_modeling.md

@@ -9,7 +9,7 @@ using ModelingToolkit, OrdinaryDiffEq
 
 @parameters t σ ρ β
 @variables x(t) y(t) z(t)
-@derivatives D'~t
+D = Differential(t)
 
 eqs = [D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
@@ -58,7 +58,7 @@ using ModelingToolkit
 
 @parameters t σ ρ β
 @variables x(t) y(t) z(t)
-@derivatives D'~t
+D = Differential(t)
 ```
 
 Then we build the system:

docs/src/tutorials/symbolic_functions.md

@@ -282,11 +282,11 @@ sufficiently large!)
 One common thing to compute in a symbolic system is derivatives. In
 ModelingToolkit.jl, derivatives are represented lazily via operations,
 just like any other function. To build a differential operator, use
-`@derivatives` like:
+`Differential` like:
 
 ```julia
 @variables t
-@derivatives D'~t
+Differential(t)
 ```
 
 This is the differential operator ``D = \frac{\partial}{\partial t}``: the number of
@@ -450,7 +450,7 @@ z = g(x) + g(y)
 
 One of the benefits of a one-language Julia symbolic stack is that the
 primitives are all written in Julia, and therefore it's trivially
-extendable from Julia itself. By default, new functions are traced
+extendible from Julia itself. By default, new functions are traced
 to the primitives and the symbolic expressions are written on the
 primitives. However, we can expand the allowed primitives by registering
 new functions. For example, let's register a new function `h`:
@@ -479,7 +479,8 @@ ModelingToolkit.derivative(::typeof(h), args::NTuple{2,Any}, ::Val{2}) = 1
 and now it works with the rest of the system:
 
 ```julia
-@derivatives Dx'~x Dy'~y
+Dx = Differential(x)
+Dy = Differential(y)
 expand_derivatives(Dx(h(x,y) + y^2)) # 2x
 expand_derivatives(Dy(h(x,y) + y^2)) # 1 + 2y
 ```

src/differentials.jl

@@ -18,8 +18,14 @@ julia> @variables x y;
 julia> D = Differential(x)
 (D'~x)
 
-julia> D(y)  # Differentiate y wrt. x
+julia> D(y) # Differentiate y wrt. x
 (D'~x)(y)
+
+julia> Dx = Differential(x) * Differential(y) # d^2/dxy operator
+(D'~x(t)) ∘ (D'~y(t))
+
+julia> D3 = Differential(x)^3 # 3rd order differential operator
+(D'~x(t)) ∘ (D'~x(t)) ∘ (D'~x(t))
 ```
 """
 struct Differential <: Function

src/systems/diffeqs/odesystem.jl

@@ -13,7 +13,7 @@ using ModelingToolkit
 
 @parameters t σ ρ β
 @variables x(t) y(t) z(t)
-@derivatives D'~t
+D = Differential(t)
 
 eqs = [D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,