diff --git a/docs/Project.toml b/docs/Project.toml
index f797fafb98..22f940c882 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,6 +1,20 @@
 [deps]
+DiffEqFlux = "aae7a2af-3d4f-5e19-a356-7da93b79d9d0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
+Integrals = "de52edbc-65ea-441a-8357-d3a637375a31"
+IntegralsCubature = "c31f79ba-6e32-46d4-a52f-182a8ac42a54"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationFlux = "253f991c-a7b2-45f8-8852-8b9a9df78a86"
+OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
+OptimizationPolyalgorithms = "500b13db-7e66-49ce-bda4-eed966be6282"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 Documenter = "0.27"
diff --git a/docs/make.jl b/docs/make.jl
index 01c56501cd..abf99459eb 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,5 +1,8 @@
 using Documenter, NeuralPDE
 
+ENV["GKSwstype"] = "100"
+using Plots
+
 include("pages.jl")
 
 makedocs(sitename = "NeuralPDE.jl",
diff --git a/docs/pages.jl b/docs/pages.jl
index d8910f0674..03d5ad35f4 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -8,10 +8,9 @@ pages = [
                                                        "pinn/low_level.md",
                                                        "pinn/ks.md",
                                                        "pinn/fp.md",
-                                                       "pinn/parm_estim.md",
+                                                       "pinn/param_estim.md",
                                                        "pinn/heterogeneous.md",
                                                        "pinn/integro_diff.md",
-                                                       "pinn/debugging.md",
                                                        "pinn/neural_adapter.md"],
     "Specialized Neural PDE Tutorials" => Any["examples/kolmogorovbackwards.md",
                                               "examples/optimal_stopping_american.md"],
@@ -23,4 +22,5 @@ pages = [
                                "solvers/kolmogorovbackwards_solver.md",
                                "solvers/optimal_stopping.md",#TODO
                                "solvers/nnrode.md"],
+    "Developer Documentation" => Any["developer/debugging.md"],
 ]
diff --git a/docs/src/pinn/debugging.md b/docs/src/developer/debugging.md
similarity index 99%
rename from docs/src/pinn/debugging.md
rename to docs/src/developer/debugging.md
index 61ef28b97c..db220cb2be 100644
--- a/docs/src/pinn/debugging.md
+++ b/docs/src/developer/debugging.md
@@ -1,5 +1,7 @@
 # Debugging PINN Solutions
 
+#### Note this is all not current right now!
+
 Let's walk through debugging functions for the physics-informed neural network
 PDE solvers.
 
diff --git a/docs/src/pinn/2D.md b/docs/src/pinn/2D.md
index be3de6b49e..0eb9613ce1 100644
--- a/docs/src/pinn/2D.md
+++ b/docs/src/pinn/2D.md
@@ -27,9 +27,8 @@ x \in [0, 2] \, ,\ y \in [0, 2] \, , \ t \in [0, 2] \, ,
 with physics-informed neural networks.
 
 ```julia
-using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
-using Quadrature, Cuba, CUDA, QuasiMonteCarlo
-import ModelingToolkit: Interval, infimum, supremum
+using NeuralPDE, Flux, Optimization, OptimizationOptimJL, DiffEqFlux
+import ModelingToolkit: Interval
 
 @parameters t x y
 @variables u(..)
diff --git a/docs/src/pinn/3rd.md b/docs/src/pinn/3rd.md
index 50cfedc0e1..a305f9aa19 100644
--- a/docs/src/pinn/3rd.md
+++ b/docs/src/pinn/3rd.md
@@ -14,7 +14,7 @@ x &\in [0, 1] \, ,
 
 We will use physics-informed neural networks.
 
-```julia
+```@example 3rdDerivative
 using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
 import ModelingToolkit: Interval, infimum, supremum
 
@@ -52,7 +52,7 @@ phi = discretization.phi
 
 We can plot the predicted solution of the ODE and its analytical solution.
 
-```julia
+```@example 3rdDerivative
 using Plots
 
 analytic_sol_func(x) = (π*x*(-x+(π^2)*(2*x-3)+1)-sin(π*x))/(π^3)
diff --git a/docs/src/pinn/fp.md b/docs/src/pinn/fp.md
index 276de0d2af..7f99f949da 100644
--- a/docs/src/pinn/fp.md
+++ b/docs/src/pinn/fp.md
@@ -20,8 +20,9 @@ p(-2.2) = p(2.2) = 0
 
 with Physics-Informed Neural Networks.
 
-```julia
+```@example fokkerplank
 using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
+using Integrals, IntegralsCubature
 import ModelingToolkit: Interval, infimum, supremum
 # the example is taken from this article https://arxiv.org/abs/1910.10503
 @parameters x
@@ -69,13 +70,15 @@ discretization = PhysicsInformedNN(chain,
                                    init_params = initθ,
                                    additional_loss=norm_loss_function)
 
-@named pde_system = PDESystem(eq,bcs,domains,[x],[p(x)])
-prob = discretize(pde_system,discretization)
+@named pdesystem = PDESystem(eq,bcs,domains,[x],[p(x)])
+prob = discretize(pdesystem,discretization)
+phi = discretization.phi
 
-pde_inner_loss_functions = prob.f.f.loss_function.pde_loss_function.pde_loss_functions.contents
-bcs_inner_loss_functions = prob.f.f.loss_function.bcs_loss_function.bc_loss_functions.contents
+sym_prob = NeuralPDE.symbolic_discretize(pdesystem, discretization)
 
-phi = discretization.phi
+pde_inner_loss_functions = sym_prob.loss_functions.pde_loss_functions
+bcs_inner_loss_functions = sym_prob.loss_functions.bc_loss_functions
+aprox_derivative_loss_functions = sym_prob.loss_functions.bc_loss_functions
 
 cb_ = function (p,l)
     println("loss: ", l )
@@ -92,7 +95,7 @@ res = Optimization.solve(prob,BFGS(),callback = cb_,maxiters=2000)
 
 And some analysis:
 
-```julia
+```@example fokkerplank
 using Plots
 C = 142.88418699042 #fitting param
 analytic_sol_func(x) = C*exp((1/(2*_σ^2))*(2*α*x^2 - β*x^4))
diff --git a/docs/src/pinn/heterogeneous.md b/docs/src/pinn/heterogeneous.md
index 59e252c041..b792d94548 100644
--- a/docs/src/pinn/heterogeneous.md
+++ b/docs/src/pinn/heterogeneous.md
@@ -1,6 +1,6 @@
-# Differential Equations with Heterogeneous Inputs
+# Differential Equations with Heterogeneous Domains
 
-A differential equation is said to have heterogeneous inputs when its dependent variables depend on different independent variables:
+A differential equation is said to have heterogeneous domains when its dependent variables depend on different independent variables:
 
 ```math
 u(x) + w(x, v) = \frac{\partial w(x, v)}{\partial w}
@@ -8,7 +8,10 @@ u(x) + w(x, v) = \frac{\partial w(x, v)}{\partial w}
 
 Here, we write an arbitrary heterogeneous system:
 
-```julia
+```@example heterogeneous
+using NeuralPDE, DiffEqFlux, Flux, ModelingToolkit, Optimization, OptimizationOptimJL
+import ModelingToolkit: Interval
+
 @parameters x y
 @variables p(..) q(..) r(..) s(..)
 Dx = Differential(x)
@@ -29,9 +32,15 @@ domains = [x ∈ Interval(0.0, 1.0),
 numhid = 3
 fastchains = [[FastChain(FastDense(1, numhid, Flux.σ), FastDense(numhid, numhid, Flux.σ), FastDense(numhid, 1)) for i in 1:2];
                         [FastChain(FastDense(2, numhid, Flux.σ), FastDense(numhid, numhid, Flux.σ), FastDense(numhid, 1)) for i in 1:2]]
-discretization = NeuralPDE.PhysicsInformedNN(fastchains, QuadratureTraining()
+discretization = NeuralPDE.PhysicsInformedNN(fastchains, QuadratureTraining())
 
 @named pde_system = PDESystem(eq, bcs, domains, [x,y], [p(x), q(y), r(x, y), s(y, x)])
 prob = SciMLBase.discretize(pde_system, discretization)
+
+callback = function (p,l)
+    println("Current loss is: $l")
+    return false
+end
+
 res = Optimization.solve(prob, BFGS(); callback = callback, maxiters=100)
 ```
diff --git a/docs/src/pinn/integro_diff.md b/docs/src/pinn/integro_diff.md
index 0a2c15dba0..fe7701268e 100644
--- a/docs/src/pinn/integro_diff.md
+++ b/docs/src/pinn/integro_diff.md
@@ -44,7 +44,7 @@ and boundary condition
 u(0) = 0
 ```
 
-```julia
+```@example integro
 using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux, DomainSets
 import ModelingToolkit: Interval, infimum, supremum
 
@@ -75,7 +75,7 @@ res = Optimization.solve(prob, BFGS(); callback = callback, maxiters=100)
 
 Plotting the final solution and analytical solution
 
-```julia
+```@example integro
 ts = [infimum(d.domain):0.01:supremum(d.domain) for d in domains][1]
 phi = discretization.phi
 u_predict  = [first(phi([t],res.minimizer)) for t in ts]
diff --git a/docs/src/pinn/ks.md b/docs/src/pinn/ks.md
index ad656790bb..725480ddbb 100644
--- a/docs/src/pinn/ks.md
+++ b/docs/src/pinn/ks.md
@@ -26,7 +26,7 @@ where `\theta = 1 - x/2` and with initial and boundary conditions:
 
 We use physics-informed neural networks.
 
-```julia
+```@example ks
 using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
 import ModelingToolkit: Interval, infimum, supremum
 
@@ -77,7 +77,7 @@ phi = discretization.phi
 
 And some analysis:
 
-```julia
+```@example ks
 using Plots
 
 xs,ts = [infimum(d.domain):dx:supremum(d.domain) for (d,dx) in zip(domains,[dx/10,dt])]
diff --git a/docs/src/pinn/parm_estim.md b/docs/src/pinn/param_estim.md
similarity index 95%
rename from docs/src/pinn/parm_estim.md
rename to docs/src/pinn/param_estim.md
index 2205ae6d97..8ee1c8ffa4 100644
--- a/docs/src/pinn/parm_estim.md
+++ b/docs/src/pinn/param_estim.md
@@ -1,6 +1,6 @@
 # Optimising Parameters of a Lorenz System
 
-Consider a Lorenz System ,
+Consider a Lorenz System,
 
 ```math
 \begin{align*}
@@ -14,7 +14,7 @@ with Physics-Informed Neural Networks. Now we would consider the case where we w
 
 We start by defining the the problem,
 
-```julia
+```@example param_estim
 using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux, OrdinaryDiffEq, Plots
 import ModelingToolkit: Interval, infimum, supremum
 @parameters t ,σ_ ,β, ρ
@@ -31,7 +31,7 @@ dt = 0.01
 
 And the neural networks as,
 
-```julia
+```@example param_estim
 input_ = length(domains)
 n = 8
 chain1 = FastChain(FastDense(input_,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,1))
@@ -43,7 +43,7 @@ We will add an additional loss term based on the data that we have in order to o
 
 Here we simply calculate the solution of the lorenz system with [OrdinaryDiffEq.jl](https://diffeq.sciml.ai/v1.10/tutorials/ode_example.html#In-Place-Updates-1) based on the adaptivity of the ODE solver. This is used to introduce non-uniformity to the time series.
 
-```julia
+```@example param_estim
 function lorenz!(du,u,p,t)
  du[1] = 10.0*(u[2]-u[1])
  du[2] = u[1]*(28.0-u[3]) - u[2]
@@ -71,7 +71,7 @@ len = length(data[2])
 
 Then we define the additional loss funciton `additional_loss(phi, θ , p)`, the function has three arguments, `phi` the trial solution, `θ` the parameters of neural networks, and the hyperparameters `p` .
 
-```julia
+```@example param_estim
 function additional_loss(phi, θ , p)
     return sum(sum(abs2, phi[i](t_ , θ[sep[i]]) .- u_[[i], :])/len for i in 1:1:3)
 end
@@ -79,7 +79,7 @@ end
 
 Then finally defining and optimising using the `PhysicsInformedNN` interface.
 
-```julia
+```@example param_estim
 discretization = NeuralPDE.PhysicsInformedNN([chain1 , chain2, chain3],NeuralPDE.GridTraining(dt), param_estim=true, additional_loss=additional_loss)
 @named pde_system = PDESystem(eqs,bcs,domains,[t],[x(t), y(t), z(t)],[σ_, ρ, β], defaults=Dict([p .=> 1.0 for p in [σ_, ρ, β]]))
 prob = NeuralPDE.discretize(pde_system,discretization)
@@ -93,7 +93,7 @@ p_ = res.minimizer[end-2:end] # p_ = [9.93, 28.002, 2.667]
 
 And then finally some analyisis by plotting.
 
-```julia
+```@example param_estim
 initθ = discretization.init_params
 acum =  [0;accumulate(+, length.(initθ))]
 sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]
diff --git a/docs/src/pinn/poisson.md b/docs/src/pinn/poisson.md
index 10b3b8e953..74e1d8d3ec 100644
--- a/docs/src/pinn/poisson.md
+++ b/docs/src/pinn/poisson.md
@@ -27,9 +27,9 @@ with grid discretization `dx = 0.1` using physics-informed neural networks.
 
 ## Copy-Pastable Code
 
-```julia
-using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
-import ModelingToolkit: Interval, infimum, supremum
+```@example
+using NeuralPDE, Flux, Optimization, OptimizationOptimJL, DiffEqFlux
+import ModelingToolkit: Interval
 
 @parameters x y
 @variables u(..)
@@ -90,9 +90,9 @@ plot(p1,p2,p3)
 
 The ModelingToolkit PDE interface for this example looks like this:
 
-```julia
+```@example poisson
 using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
-import ModelingToolkit: Interval, infimum, supremum
+import ModelingToolkit: Interval
 
 @parameters x y
 @variables u(..)
@@ -113,7 +113,7 @@ domains = [x ∈ Interval(0.0,1.0),
 Here, we define the neural network, where the input of NN equals the number of dimensions and output equals the number of equations in the system.
 
 
-```julia
+```@example poisson
 # Neural network
 dim = 2 # number of dimensions
 chain = FastChain(FastDense(dim,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))
@@ -121,7 +121,7 @@ chain = FastChain(FastDense(dim,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(1
 
 Convert weights of neural network from Float32 to Float64 in order to all inner calculation will be with Float64.
 
-```julia
+```@example poisson
 # Initial parameters of Neural network
 initθ = Float64.(DiffEqFlux.initial_params(chain))
 ```
@@ -129,7 +129,7 @@ initθ = Float64.(DiffEqFlux.initial_params(chain))
 Here, we build PhysicsInformedNN algorithm where `dx` is the step of discretization, `strategy` stores information for choosing a training strategy and
 `init_params =initθ` initial parameters of neural network.
 
-```julia
+```@example poisson
 # Discretization
 dx = 0.05
 discretization = PhysicsInformedNN(chain, GridTraining(dx),init_params =initθ)
@@ -137,7 +137,7 @@ discretization = PhysicsInformedNN(chain, GridTraining(dx),init_params =initθ)
 
 As described in the API docs, we now need to define the `PDESystem` and create PINNs problem using the `discretize` method.
 
-```julia
+```@example poisson
 @named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
 prob = discretize(pde_system,discretization)
 ```
@@ -145,7 +145,7 @@ prob = discretize(pde_system,discretization)
 Here, we define the callback function and the optimizer. And now we can solve the PDE using PINNs
 (with the number of epochs `maxiters=1000`).
 
-```julia
+```@example poisson
 #Optimizer
 opt = OptimizationOptimJL.BFGS()
 
@@ -160,7 +160,7 @@ phi = discretization.phi
 
 We can plot the predicted solution of the PDE and compare it with the analytical solution in order to plot the relative error.
 
-```julia
+```@example poisson
 xs,ys = [infimum(d.domain):dx/10:supremum(d.domain) for d in domains]
 analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)
 
@@ -168,6 +168,8 @@ u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(l
 u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
 diff_u = abs.(u_predict .- u_real)
 
+using Plots
+
 p1 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
 p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
 p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
diff --git a/docs/src/pinn/system.md b/docs/src/pinn/system.md
index 24076c653e..3730e4ca7f 100644
--- a/docs/src/pinn/system.md
+++ b/docs/src/pinn/system.md
@@ -32,9 +32,10 @@ u_2(t, 0) & = - u_2(t, 1) = e^{-t} \, ,
 
 with physics-informed neural networks.
 
-```julia
+## Solution
+
+```@example system
 using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
-using Quadrature,Cubature
 import ModelingToolkit: Interval, infimum, supremum
 
 @parameters t, x
@@ -68,18 +69,18 @@ n = 15
 chain =[FastChain(FastDense(input_,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,1)) for _ in 1:3]
 initθ = map(c -> Float64.(c), DiffEqFlux.initial_params.(chain))
 
-_strategy = QuadratureTraining()
-discretization = PhysicsInformedNN(chain, _strategy, init_params= initθ)
+strategy = QuadratureTraining()
+discretization = PhysicsInformedNN(chain, strategy, init_params = initθ)
 
-@named pde_system = PDESystem(eqs,bcs,domains,[t,x],[u1(t, x),u2(t, x),u3(t, x)])
-prob = discretize(pde_system,discretization)
-sym_prob = symbolic_discretize(pde_system,discretization)
+@named pdesystem = PDESystem(eqs,bcs,domains,[t,x],[u1(t, x),u2(t, x),u3(t, x)])
+prob = discretize(pdesystem,discretization)
+sym_prob = symbolic_discretize(pdesystem,discretization)
 
-pde_inner_loss_functions = prob.f.f.loss_function.pde_loss_function.pde_loss_functions.contents
-bcs_inner_loss_functions = prob.f.f.loss_function.bcs_loss_function.bc_loss_functions.contents
+pde_inner_loss_functions = sym_prob.loss_functions.pde_loss_functions
+bcs_inner_loss_functions = sym_prob.loss_functions.bc_loss_functions
 
-callback = function (p,l)
-    println("loss: ", l )
+callback = function (p, l)
+    println("loss: ", l)
     println("pde_losses: ", map(l_ -> l_(p), pde_inner_loss_functions))
     println("bcs_losses: ", map(l_ -> l_(p), bcs_inner_loss_functions))
     return false
@@ -90,11 +91,14 @@ res = Optimization.solve(prob,BFGS(); callback = callback, maxiters=5000)
 phi = discretization.phi
 ```
 
-Low-level api
+## Direct Construction via symbolic_discretize
+
+One can take apart the pieces and reassemble the loss functions using the `symbolic_discretize`
+interface. Here is an example using the components from `symbolic_discretize` to fully
+reproduce the `discretize` optimization:
 
-```julia
+```@example system
 using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
-using Quadrature,Cubature
 import ModelingToolkit: Interval, infimum, supremum
 
 @parameters t, x
@@ -129,51 +133,21 @@ initθ = map(c -> Float64.(c), DiffEqFlux.initial_params.(chain))
 flat_initθ = reduce(vcat,initθ )
 
 eltypeθ = eltype(initθ[1])
-phi = NeuralPDE.get_phi.(chain)
-
-map(phi_ -> phi_(rand(2,10), flat_initθ),phi)
-
-derivative = NeuralPDE.get_numeric_derivative()
-
-
-indvars = [t,x]
-depvars = [u1,u2,u3]
-dim = length(domains)
-quadrature_strategy = NeuralPDE.QuadratureTraining()
-
+@named pdesystem = PDESystem(eqs,bcs,domains,[t,x],[u1(t, x),u2(t, x),u3(t, x)])
 
-_pde_loss_functions = [NeuralPDE.build_loss_function(eq,indvars,depvars,phi,derivative,
-                                                     chain,initθ,quadrature_strategy) for eq in  eqs]
+strategy = NeuralPDE.QuadratureTraining()
+discretization = PhysicsInformedNN(chain, strategy, init_params = initθ)
+sym_prob = NeuralPDE.symbolic_discretize(pdesystem, discretization)
 
-map(loss_f -> loss_f(rand(2,10), flat_initθ),_pde_loss_functions)
+pde_loss_functions = sym_prob.loss_functions.pde_loss_functions
+bc_loss_functions = sym_prob.loss_functions.bc_loss_functions
 
-bc_indvars = NeuralPDE.get_argument(bcs,indvars,depvars)
-_bc_loss_functions = [NeuralPDE.build_loss_function(bc,indvars,depvars, phi, derivative,
-                                                    chain,initθ,quadrature_strategy,
-                                                    bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]
-map(loss_f -> loss_f(rand(1,10), flat_initθ),_bc_loss_functions)
-
-# dx = 0.1
-# train_sets = NeuralPDE.generate_training_sets(domains,dx,eqs,bcs,eltypeθ,indvars,depvars)
-# pde_train_set,bcs_train_set = train_sets
-pde_bounds, bcs_bounds = NeuralPDE.get_bounds(domains,eqs,bcs,eltypeθ,indvars,depvars,quadrature_strategy)
-
-plbs,pubs = pde_bounds
-pde_loss_functions = [NeuralPDE.get_loss_function(_loss,
-                                                 lb,ub,
-                                                 eltypeθ,
-                                                 quadrature_strategy)
-                                                 for (_loss,lb,ub) in zip(_pde_loss_functions, plbs,pubs)]
-
-map(l->l(flat_initθ) ,pde_loss_functions)
-
-blbs,bubs = bcs_bounds
-bc_loss_functions = [NeuralPDE.get_loss_function(_loss,lb,ub,
-                                                 eltypeθ,
-                                                 quadrature_strategy)
-                                                 for (_loss,lb,ub) in zip(_bc_loss_functions, blbs,bubs)]
-
-map(l->l(flat_initθ) ,bc_loss_functions)
+callback = function (p, l)
+    println("loss: ", l)
+    println("pde_losses: ", map(l_ -> l_(p), pde_loss_functions))
+    println("bcs_losses: ", map(l_ -> l_(p), bc_loss_functions))
+    return false
+end
 
 loss_functions =  [pde_loss_functions;bc_loss_functions]
 
@@ -184,21 +158,15 @@ end
 f_ = OptimizationFunction(loss_function, Optimization.AutoZygote())
 prob = Optimization.OptimizationProblem(f_, flat_initθ)
 
-cb_ = function (p,l)
-    println("loss: ", l )
-    println("pde losses: ", map(l -> l(p), loss_functions[1:3]))
-    println("bcs losses: ", map(l -> l(p), loss_functions[4:end]))
-    return false
-end
-
-res = Optimization.solve(prob,OptimizationOptimJL.BFGS(); callback = cb_, maxiters=5000)
+res = Optimization.solve(prob,OptimizationOptimJL.BFGS(); callback = callback, maxiters=5000)
 ```
 
-And some analysis for both low and high level api:
+And some analysis for both the `symbolic_discretize` and `discretize` APIs:
 
-```julia
+```@example system
 using Plots
 
+phi = discretization.phi
 ts,xs = [infimum(d.domain):0.01:supremum(d.domain) for d in domains]
 
 acum =  [0;accumulate(+, length.(initθ))]
@@ -227,16 +195,25 @@ end
 ## Derivative neural network approximation
 
 The accuracy and stability of numerical derivative decreases with each successive order.
-The accuracy of the entire solution is determined by the worst accuracy of one of the variables, in our case - the highest degree of the derivative.
-Derivative neural network approximation is such an approach that using lower-order numeric derivatives and estimates higher-order derivatives with a neural network so that allows an increase in the marginal precision for all optimization.
-
-Since `u3` is only in the first and second equations, that its accuracy during training is determined by the accuracy of the second numerical derivative `u3(t,x) ~ (Dtt(u1(t,x)) -Dxx(u1(t,x))) / sin(pi*x)`.
-
-We approximate the derivative of the neural network with another neural network `Dt(u1(t,x)) ~ Dtu1(t,x)` and train it along with other equations, and thus we avoid using the second numeric derivative `Dt(Dtu1(t,x))`.
-
-```julia
-using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
-using Quadrature,Cubature
+The accuracy of the entire solution is determined by the worst accuracy of one of the variables,
+in our case - the highest degree of the derivative. Meanwhile the computational cost of automatic
+differentation for higher orders grows the power in O(n^d), making even numerical differentiation
+much more efficient! Given these two bad choices, there exists an alternative which can improve
+training speed and accuracy: using a system to represent the derivatives directly.
+
+The derivative neural network approximation is such an approach that using lower-order numeric
+derivatives and estimates higher-order derivatives with a neural network so that allows an
+increase in the marginal precision for all optimization. Since `u3` is only in the first and
+second equations, that its accuracy during training is determined by the accuracy of the
+second numerical derivative `u3(t,x) ~ (Dtt(u1(t,x)) -Dxx(u1(t,x))) / sin(pi*x)`.
+
+We approximate the derivative of the neural network with another neural network
+`Dt(u1(t,x)) ~ Dtu1(t,x)` and train it along with other equations, and thus we avoid
+using the second numeric derivative `Dt(Dtu1(t,x))`.
+
+```@example derivativenn
+using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationFlux, OptimizationOptimJL, DiffEqFlux
+using Plots
 import ModelingToolkit: Interval, infimum, supremum
 
 @parameters t, x
@@ -267,6 +244,10 @@ der_ = [Dt(u1(t,x)) ~ Dtu1(t,x),
 
 bcs__ = [bcs_;der_]
 
+# Space and time domains
+domains = [t ∈ Interval(0.0, 1.0),
+           x ∈ Interval(0.0, 1.0)]
+
 input_ = length(domains)
 n = 15
 chain = [FastChain(FastDense(input_,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,1)) for _ in 1:7]
@@ -278,14 +259,13 @@ discretization = NeuralPDE.PhysicsInformedNN(chain,
                                              init_params= initθ)
 
 vars = [u1(t,x),u2(t,x),u3(t,x),Dxu1(t,x),Dtu1(t,x),Dxu2(t,x),Dtu2(t,x)]
-@named pde_system = PDESystem(eqs_,bcs__,domains,[t,x],vars)
-prob = NeuralPDE.discretize(pde_system,discretization)
-sym_prob = NeuralPDE.symbolic_discretize(pde_system,discretization)
+@named pdesystem = PDESystem(eqs_,bcs__,domains,[t,x],vars)
+prob = NeuralPDE.discretize(pdesystem,discretization)
+sym_prob = NeuralPDE.symbolic_discretize(pdesystem,discretization)
 
-pde_inner_loss_functions = prob.f.f.loss_function.pde_loss_function.pde_loss_functions.contents
-inner_loss_functions = prob.f.f.loss_function.bcs_loss_function.bc_loss_functions.contents
-bcs_inner_loss_functions = inner_loss_functions[1:7]
-aprox_derivative_loss_functions = inner_loss_functions[9:end]
+pde_inner_loss_functions = sym_prob.loss_functions.pde_loss_functions
+bcs_inner_loss_functions = sym_prob.loss_functions.bc_loss_functions[1:7]
+aprox_derivative_loss_functions = sym_prob.loss_functions.bc_loss_functions[9:end]
 
 callback = function (p,l)
     println("loss: ", l )
@@ -304,12 +284,11 @@ phi = discretization.phi
 
 And some analysis:
 
-```julia
+```@example derivativenn
 using Plots
 
 ts,xs = [infimum(d.domain):0.01:supremum(d.domain) for d in domains]
 
-initθ = discretization.init_params
 acum =  [0;accumulate(+, length.(initθ))]
 sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]
 minimizers_ = [res.minimizer[s] for s in sep]
@@ -354,17 +333,20 @@ end
 
 Also, in addition to systems, we can use the matrix form of PDEs:
 
-```julia
+```@example
+using ModelingToolkit, NeuralPDE
 @parameters x y
 @variables u[1:2,1:2](..)
 @derivatives Dxx''~x
 @derivatives Dyy''~y
 
+# Initial and boundary conditions
+bcs = [u[1](x,0) ~ x, u[2](x,0) ~ 2, u[3](x,0) ~ 3, u[4](x,0) ~ 4]
+
 # matrix PDE
 eqs  = @. [(Dxx(u_(x,y)) + Dyy(u_(x,y))) for u_ in u] ~ -sin(pi*x)*sin(pi*y)*[0 1; 0 1]
 
-# Initial and boundary conditions
-bcs = [u[1](x,0) ~ x, u[2](x,0) ~ 2, u[3](x,0) ~ 3, u[4](x,0) ~ 4]
+size(eqs)
 ```
 
 ## Linear parabolic system of PDEs
@@ -392,10 +374,9 @@ w(t, 1) = \frac{e^{\lambda_1} cos(\frac{x}{a})-e^{\lambda_2}cos(\frac{x}{a})}{\l
 
 with a physics-informed neural network.
 
-```julia
+```@example
 using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
 using Plots
-using Quadrature,Cubature
 import ModelingToolkit: Interval, infimum, supremum
 
 @parameters t, x
@@ -439,15 +420,15 @@ n = 15
 chain = [FastChain(FastDense(input_, n, Flux.σ), FastDense(n, n, Flux.σ), FastDense(n, 1)) for _ in 1:2]
 initθ = map(c -> Float64.(c), DiffEqFlux.initial_params.(chain))
 
-_strategy = QuadratureTraining()
-discretization = PhysicsInformedNN(chain, _strategy, init_params=initθ)
+strategy = QuadratureTraining()
+discretization = PhysicsInformedNN(chain, strategy, init_params=initθ)
 
-@named pde_system = PDESystem(eqs, bcs, domains, [t,x], [u(t,x),w(t,x)])
-prob = discretize(pde_system, discretization)
-sym_prob = symbolic_discretize(pde_system, discretization)
+@named pdesystem = PDESystem(eqs, bcs, domains, [t,x], [u(t,x),w(t,x)])
+prob = discretize(pdesystem, discretization)
+sym_prob = symbolic_discretize(pdesystem, discretization)
 
-pde_inner_loss_functions = prob.f.f.loss_function.pde_loss_function.pde_loss_functions.contents
-bcs_inner_loss_functions = prob.f.f.loss_function.bcs_loss_function.bc_loss_functions.contents
+pde_inner_loss_functions = sym_prob.loss_functions.pde_loss_functions
+bcs_inner_loss_functions = sym_prob.loss_functions.bc_loss_functions
 
 callback = function (p, l)
     println("loss: ", l)
@@ -511,10 +492,9 @@ where k is a root of the algebraic (transcendental) equation f(k) = g(k).
 
 This is done using a derivative neural network approximation.
 
-```julia
-using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux, DifferentialEquations, Roots
+```@example
+using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux, Roots
 using Plots
-using Quadrature,Cubature
 import ModelingToolkit: Interval, infimum, supremum
 
 @parameters x, y
@@ -565,18 +545,21 @@ n = 15
 chain = [FastChain(FastDense(input_, n, Flux.σ), FastDense(n, n, Flux.σ), FastDense(n, 1)) for _ in 1:6] # 1:number of @variables
 initθ = map(c -> Float64.(c), DiffEqFlux.initial_params.(chain))
 
-_strategy = QuadratureTraining()
-discretization = PhysicsInformedNN(chain, _strategy, init_params=initθ)
+strategy = QuadratureTraining()
+discretization = PhysicsInformedNN(chain, strategy, init_params=initθ)
 
 vars = [u(x,y),w(x,y),Dxu(x,y),Dyu(x,y),Dxw(x,y),Dyw(x,y)]
-@named pde_system = PDESystem(eqs_, bcs__, domains, [x,y], vars)
-prob = NeuralPDE.discretize(pde_system, discretization)
-sym_prob = NeuralPDE.symbolic_discretize(pde_system, discretization)
+@named pdesystem  = PDESystem(eqs_, bcs__, domains, [x,y], vars)
+prob = NeuralPDE.discretize(pdesystem , discretization)
+sym_prob = NeuralPDE.symbolic_discretize(pdesystem , discretization)
+
+strategy = NeuralPDE.QuadratureTraining()
+discretization = PhysicsInformedNN(chain, strategy, init_params = initθ)
+sym_prob = NeuralPDE.symbolic_discretize(pdesystem, discretization)
 
-pde_inner_loss_functions = prob.f.f.loss_function.pde_loss_function.pde_loss_functions.contents
-inner_loss_functions = prob.f.f.loss_function.bcs_loss_function.bc_loss_functions.contents
-bcs_inner_loss_functions = inner_loss_functions[1:6]
-aprox_derivative_loss_functions = inner_loss_functions[7:end]
+pde_inner_loss_functions = sym_prob.loss_functions.pde_loss_functions
+bcs_inner_loss_functions = sym_prob.loss_functions.bc_loss_functions[1:6]
+aprox_derivative_loss_functions = sym_prob.loss_functions.bc_loss_functions[7:end]
 
 callback = function (p, l)
     println("loss: ", l)
@@ -647,11 +630,10 @@ where k is a root of the algebraic (transcendental) equation f(k) = g(k), j0 and
 
 We solve this with Neural:
 
-```julia
+```@example
 using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux, Roots
 using SpecialFunctions
 using Plots
-using Quadrature,Cubature
 import ModelingToolkit: Interval, infimum, supremum
 
 @parameters t, x
@@ -671,7 +653,7 @@ g(x) = 4 * cos(π * x)
 root(x) = g(x) - f(x)
 
 # Analytic solution
-k = find_zero(root, (0, 1), Bisection())                # k is a root of the algebraic (transcendental) equation f(x) = g(x)
+k = find_zero(root, (0, 1), Roots.Bisection())                # k is a root of the algebraic (transcendental) equation f(x) = g(x)
 ξ(t, x) = sqrt(f(k)) / sqrt(a) * sqrt(a * (t + 1)^2 - (x + 1)^2)
 θ(t, x) = besselj0(ξ(t, x)) + bessely0(ξ(t, x))                     # Analytical solution to Klein-Gordon equation
 w_analytic(t, x) = θ(t, x)
@@ -699,15 +681,15 @@ n = 15
 chain = [FastChain(FastDense(input_, n, Flux.σ), FastDense(n, n, Flux.σ), FastDense(n, 1)) for _ in 1:2]
 initθ = map(c -> Float64.(c), DiffEqFlux.initial_params.(chain))
 
-_strategy = QuadratureTraining()
-discretization = PhysicsInformedNN(chain, _strategy, init_params=initθ)
+strategy = QuadratureTraining()
+discretization = PhysicsInformedNN(chain, strategy, init_params=initθ)
 
-@named pde_system = PDESystem(eqs, bcs, domains, [t,x], [u(t,x),w(t,x)])
-prob = discretize(pde_system, discretization)
-sym_prob = symbolic_discretize(pde_system, discretization)
+@named pdesystem = PDESystem(eqs, bcs, domains, [t,x], [u(t,x),w(t,x)])
+prob = discretize(pdesystem, discretization)
+sym_prob = symbolic_discretize(pdesystem, discretization)
 
-pde_inner_loss_functions = prob.f.f.loss_function.pde_loss_function.pde_loss_functions.contents
-bcs_inner_loss_functions = prob.f.f.loss_function.bcs_loss_function.bc_loss_functions.contents
+pde_inner_loss_functions = sym_prob.loss_functions.pde_loss_functions
+bcs_inner_loss_functions = sym_prob.loss_functions.bc_loss_functions
 
 callback = function (p, l)
     println("loss: ", l)
diff --git a/docs/src/pinn/wave.md b/docs/src/pinn/wave.md
index 182b84bdba..a4cfa08ceb 100644
--- a/docs/src/pinn/wave.md
+++ b/docs/src/pinn/wave.md
@@ -15,9 +15,9 @@ with grid discretization `dx = 0.1` and physics-informed neural networks.
 
 Further, the solution of this equation with the given boundary conditions is presented.
 
-```julia
-using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
-import ModelingToolkit: Interval, infimum, supremum
+```@example wave
+using NeuralPDE, Flux, Optimization, OptimizationOptimJL, DiffEqFlux
+import ModelingToolkit: Interval
 
 @parameters t, x
 @variables u(..)
@@ -62,7 +62,7 @@ phi = discretization.phi
 
 We can plot the predicted solution of the PDE and compare it with the analytical solution in order to plot the relative error.
 
-```julia
+```@example wave
 using Plots
 
 ts,xs = [infimum(d.domain):dx:supremum(d.domain) for d in domains]
@@ -95,7 +95,7 @@ u_t(0, x) = 1 - 2x \\
 
 with grid discretization `dx = 0.05` and physics-informed neural networks. Here we take advantage of adaptive derivative to increase accuracy.
 
-```julia
+```@example wave2
 using NeuralPDE, Flux, ModelingToolkit, Optimization, OptimizationOptimJL, DiffEqFlux
 using Plots, Printf
 import ModelingToolkit: Interval, infimum, supremum
@@ -151,13 +151,13 @@ discretization = PhysicsInformedNN(chain, strategy; init_params=initθ)
 
 @named pde_system = PDESystem(eq, bcs, domains, [t, x], [u(t, x), Dxu(t, x), Dtu(t, x), O1(t, x), O2(t, x)])
 prob = discretize(pde_system, discretization)
+sym_prob = NeuralPDE.symbolic_discretize(pde_system, discretization)
 
-pde_inner_loss_functions = prob.f.f.loss_function.pde_loss_function.pde_loss_functions.contents
-inner_loss_functions = prob.f.f.loss_function.bcs_loss_function.bc_loss_functions.contents
-bcs_inner_loss_functions = inner_loss_functions
+pde_inner_loss_functions = sym_prob.loss_functions.pde_loss_functions
+bcs_inner_loss_functions = sym_prob.loss_functions.bc_loss_functions
 
 callback = function (p, l)
-    println("Current loss is: $l")
+    println("loss: ", l)
     println("pde_losses: ", map(l_ -> l_(p), pde_inner_loss_functions))
     println("bcs_losses: ", map(l_ -> l_(p), bcs_inner_loss_functions))
     return false