Update README

Author: dlfivefifty

README.md

@@ -36,10 +36,9 @@ h = f + g^2
 r = roots(h)
 rp = roots(h')
 
-using Plots
-plot(h)
-scatter!(r,h.(r))
-scatter!(rp,h.(rp))
+plot(h; label="f + g^2")
+scatter!(r, h.(r); label="roots")
+scatter!(rp, h.(rp); label="extrema")
 ```
 
 <img src=https://github.com/JuliaApproximation/ApproxFun.jl/raw/master/images/extrema.png width=500 height=400>
@@ -54,15 +53,15 @@ derivative, the `norm` is small:
 
 ```julia
 f = Fun(x->exp(x), -1..1)
-norm(f-f')
+norm(f-f')  # 4.4391656415701095e-14
 ```
 
 Similarly, `cumsum` defines an indefinite integration operator:
 
 ```julia
 g = cumsum(f)
 g = g + f(-1)
-norm(f-g)
+norm(f-g) # 3.4989733283850415e-15d
 ```
 
 Algebraic and differential operations are also implemented where possible, and most of Julia's built-in functions are overridden to accept `Fun`s:
@@ -78,17 +77,15 @@ h = airyai(10asin(f)+2g)
 ## Solving ordinary differential equations
 
 
-Solve the Airy ODE `u'' - x u = 0` as a BVP on `[-1000,200]`:
+We can also solve differential equations. Consider the Airy ODE `u'' - x u = 0` as a boundary value problem on `[-1000,200]` with conditions `u(-1000) = 1` and `u(200) = 2`. The unique solution is a linear combination of Airy functions. We can calculate it as follows:
 
 ```julia
-a,b = -1000,200
-x = Fun(identity, a..b)
-d = domain(x)
-D = Derivative(d)
-B = Dirichlet(d)
-L = D^2 - x
-u = [B;L] \ [[airyai(a),airyai(b)],0]
-plot(u)
+x = Fun(identity, -1000..200) # the function x on the interval -1000..200
+D = Derivative()              # The derivative operator
+B = Dirichlet()               # Dirichlet conditions
+L = D^2 - x                   # the Airy operator
+u = [B;L] \ [[1,2],0]         # Calculate u such that B*u == [1,2] and L*u == 0
+plot(u; label="u")
 ```
 
 <img src=https://github.com/JuliaApproximation/ApproxFun.jl/raw/master/images/airy.png width=500 height=400>
@@ -100,10 +97,9 @@ Solve a nonlinear boundary value problem satisfying the ODE `0.001u'' + 6*(1-x^2
 
 ```julia
 x  = Fun()
-u₀ = 0.0x
-
+u₀ = 0.0x # initial guess
 N = u -> [u(-1)-1, u(1)+0.5, 0.001u'' + 6*(1-x^2)*u' + u^2 - 1]
-u = newton(N,u0)
+u = newton(N, u₀) # perform Newton iteration in function space
 plot(u)
 ```
 
@@ -112,22 +108,18 @@ plot(u)
 One can also solve a system nonlinear ODEs with potentially nonlinear boundary conditions:
 
 ```julia
-function nonlinear_test()
-    x=Fun(Chebyshev(0..1))
-    N = (u1,u2) -> [u1'(0) - 0.5*u1(0)*u2(0);
-                    u2'(0) + 1;
-                    u1(1) - 1;
-                    u2(1) - 1;
-                    u1'' + u1*u2;
-                    u2'' - u1*u2]
-
-    u10 = one(x)
-    u20 = one(x)
-
-    return newton(N, [u10,u20])
-end
 
-u1,u2 = nonlinear_test()
+x=Fun(identity, 0..1)
+N = (u1,u2) -> [u1'(0) - 0.5*u1(0)*u2(0);
+                u2'(0) + 1;
+                u1(1) - 1;
+                u2(1) - 1;
+                u1'' + u1*u2;
+                u2'' - u1*u2]
+
+u10 = one(x)
+u20 = one(x)
+u1,u2 = newton(N, [u10,u20])
 
 plot(u1, label="u1")
 plot!(u2, label="u2")
@@ -144,7 +136,7 @@ Specify the space `Fourier` to ensure that the representation is periodic:
 
 ```julia
 f = Fun(cos, Fourier(-π..π))
-norm(f' + Fun(sin, Fourier(-π..π))
+norm(f' + Fun(sin, Fourier(-π..π))) # 5.923502902288505e-17
 ```
 
 Due to the periodicity, Fourier representations allow for the asymptotic savings of `2/π`

images/airy.png

@@ -7,7 +7,6 @@ using ApproxFun, SpecialFunctions, LinearAlgebra, Test
         f = sin(x^2)
         g = cos(x)
 
-
         @test ≈(f(.1),sin(.1^2);atol=1000eps())
 
         h = f + g^2