more improvements

Author: ChrisRackauckas

docs/src/basics/Validation.md

@@ -83,7 +83,7 @@ second argument.
 ```@example validation2
 using ModelingToolkit, Unitful
 # Composite type parameter in registered function
-@parameters t
+@variables t
 D = Differential(t)
 struct NewType
     f::Any

docs/src/basics/Variable_metadata.md

@@ -10,14 +10,14 @@ Descriptive strings can be attached to variables using the `[description = "desc
 
 ```@example metadata
 using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
 @variables u [description = "This is my input"]
 getdescription(u)
 ```
 
 When variables with descriptions are present in systems, they will be printed when the system is shown in the terminal:
 
 ```@example metadata
-@parameters t
 @variables u(t) [description = "A short description of u"]
 @parameters p [description = "A description of p"]
 @named sys = ODESystem([u ~ p], t)
@@ -62,6 +62,8 @@ Designate a variable as either an input or an output using the following
 
 ```@example metadata
 using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
 @variables u [input = true]
 isinput(u)
 ```
@@ -137,8 +139,6 @@ For systems that contain parameters with metadata like described above, have som
 In the example below, we define a system with tunable parameters and extract bounds vectors
 
 ```@example metadata
-@parameters t
-Dₜ = Differential(t)
 @variables x(t)=0 u(t)=0 [input = true] y(t)=0 [output = true]
 @parameters T [tunable = true, bounds = (0, Inf)]
 @parameters k [tunable = true, bounds = (0, Inf)]

docs/src/examples/higher_order.md

@@ -22,7 +22,7 @@ eqs = [D(D(x)) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs,t)
 ```
 
 Note that we could've used an alternative syntax for 2nd order, i.e.

docs/src/examples/modelingtoolkitize_index_reduction.md

@@ -31,7 +31,7 @@ pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
 traced_sys = modelingtoolkitize(pendulum_prob)
 pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
 prob = ODEProblem(pendulum_sys, [], tspan)
-sol = solve(prob, Tsit5(), abstol = 1e-8, reltol = 1e-8)
+sol = solve(prob, Rodas5P(), abstol = 1e-8, reltol = 1e-8)
 plot(sol, idxs = unknowns(traced_sys))
 ```
 
@@ -71,7 +71,7 @@ u0 = [1.0, 0, 0, 0, 0];
 p = [9.8, 1];
 tspan = (0, 10.0);
 pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
-solve(pendulum_prob, Rodas4())
+solve(pendulum_prob, Rodas5P())
 ```
 
 However, one will quickly be greeted with the unfortunate message:
@@ -151,7 +151,7 @@ numerical solver. Let's try that out:
 traced_sys = modelingtoolkitize(pendulum_prob)
 pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
 prob = ODEProblem(pendulum_sys, Pair[], tspan)
-sol = solve(prob, Rodas4())
+sol = solve(prob, Rodas5P())
 
 using Plots
 plot(sol, idxs = unknowns(traced_sys))

docs/src/tutorials/ode_modeling.md

@@ -256,7 +256,7 @@ function fol_factory(separate = false; name)
         D(x) ~ RHS] :
           D(x) ~ (f - x) / τ
 
-    ODESystem(eqs; name)
+    ODESystem(eqs, t; name)
 end
 ```
 
@@ -277,7 +277,7 @@ again are just algebraic relations:
 connections = [fol_1.f ~ 1.5,
     fol_2.f ~ fol_1.x]
 
-connected = compose(ODESystem(connections, name = :connected), fol_1, fol_2)
+connected = compose(ODESystem(connections, t, name = :connected), fol_1, fol_2)
 ```
 
 All equations, variables, and parameters are collected, but the structure of the