Numerical solution to the Lane-Emden equation

Question

I am new to Julia programming, but for fun and practice, I wrote a program that numerically solves the Lane-Emden equation. I even wrote an interactive version that changes the values of n and log_delta_xi with sliders from PlutoUI (not included here).

I will leave the program here and encourage any suggestions as to how it can be improved.

using LaTeXStrings
using Plots
using Markdown
gr()

println("Lane-Emden Equation")
L"""\dfrac{\text{d}}{\text{d}\xi} \ 
    \left( \xi^2 \dfrac{\text{d}\theta}{\text{d}\xi} \right) \
    = \
    -\xi^2\theta^n \ """

display(md"""
Separation of variables

``
\dfrac{\text{d}y}{\text{d}\xi}   =  \dfrac{z}{\xi^2} 
``

``
\dfrac{\text{d}z}{\text{d}\xi}  = -\xi^2y^n
``
""")

function solveLaneEmden(log_delta_xi=-4, n=3)
    delta_xi = 10.0^log_delta_xi
    
    # Inner boundary condition 
    y0 = 1 - delta_xi^2/6 
    z0 = -delta_xi^3/3 
    
    ys  = [y0]
    zs  = [z0]
    xis = [delta_xi]
    ycs = [y0]
    zcs = [z0]
    
    while true
        y  =  last(ys)
        z  =  last(zs)
        xi =  last(xis)
        yc =  last(ycs)
        zc =  last(zcs)
        
           ## Primitive method 
        yi = y + delta_xi * z/xi^2
        zi = z + delta_xi * -xi^2*y^n
        
        ## Predictor-corrector technique 
        xii = xi + delta_xi
        yci = yc + 1/2 * delta_xi * (z/xi^2 + zi/xii^2)
        zci = zc + 1/2 * delta_xi * (-xi^2*y^n - xi^2*yi^n)
        
        # Outer boundary condition 
        if (yi < 1e-10 || yci < 1e-10)
            break
        end
        
        push!(xis, xii)
        push!(ys, yi)
        push!(zs, zi)
        push!(ycs, yci)
        push!(zcs, zci)
            
    end

    return (xis, ys, ycs)
    
end

function plotLaneEmden(log_delta_xi=-4, n=3)
    xis, ys, ycs = solveLaneEmden(log_delta_xi, n)
    
    xi2 = range(0,sqrt(6),step=1e-3)
    # @. will add . to every operator
    Plots.plot(xi2, 1 .- xi2.^2/6, label="n = 0")
    
    xi2 = range(0, pi, step=1e-3) 
    Plots.plot!(xi2, sin.(xi2)./xi2, linecolor = :orange, label="n = 1")
    
    Plots.plot!(xis, ys,  linecolor = :green, label="n = $n")
    Plots.plot!(xis, ycs, linecolor = :black, linestyle = :dash, label="P-C")
    
    Plots.xlabel!("ξ")
    Plots.ylabel!("θ")

end

plotLaneEmden()