CalculusWithJulia ODEs Part Question

General Usage
1

Hi all,

I want to know why the code in CalculusWithJulia produce different solution for the same Initial Value problem:

y' = 0.2y(16-y)
y(0)=3

The plot is great anyway, since CairoMakie is a bit difficult to handle, this one really comes in handy.

so this is my problem, why this code:

# https://jverzani.github.io/CalculusWithJuliaNotes.jl/ODEs/odes.html
using SymPy, CalculusWithJulia, LaTeXStrings,Plots
@syms x
@syms y()

println("First-Order Linear Differential Equations with Julia")
println("y' = 0.2y(16-y) is the differential equation")
println("")

#  we can use diff(u(x),x) but here, for visual simplicity, use the Differential operator,
D = Differential(x)
eqn = D(y)(x) ~ 0.2 * y(x) * (16 - y(x))
out = dsolve(eqn)
# The limits are confirmed by investigating the limits of the right-hand:
limit(rhs(out), x => oo), limit(rhs(out), x => -oo)

# For the initial value problem y(0)=3
eq = rhs(out)    # just the right hand side
C1 = first(setdiff(free_symbols(eq), (x))) # fish out constant
c1 = solve(eq(x=>0) - 3, C1)
eq(C1 => c1[1])

sol = out(C1 => c1[1])

x0, y0 = 0, 3
ivp_sol = dsolve(eqn, y(x), ics=Dict(y(x0) => y0))

println("the particular solution for y(0) = 3:  ")
sympy.pretty_print(ivp_sol)

produce different IVP solution than this code:

using SymPy, LaTeXStrings,Plots,PrettyTables
@syms x,y,c
@syms f()

println("First-Order Linear Differential Equations with Julia")
println("y' = 0.2y(16-y) is the differential equation")
println("")

p1(x) = -(3.2-0.2y)*x^0
g1(x) = 0
μ1(x) = exp(integrate(p1(x),x))
y1(x) = (1/μ1(x))*(integrate((μ1(x)*g1(x)),x) + c)

println("The general solution:  ")
sympy.pretty_print(y1(x))
println("")


# Computing IVP with SymPy
diffeqf = Eq(f(x).diff(x) -(3.2 -0.2f(x))*x^0, g1(x)); string(diffeqf)
ivp_f = dsolve(diffeqf, f(x), ics = Dict(f(0) => 3))

println("the particular solution for velocity, f(0) = 3:  ")
sympy.pretty_print(ivp_f)

I am able to plot it with this code:

using CalculusWithJulia, Plots, SymPy

@syms x y
x0, y0 = 0, 0

F(y, x) = 0.2y*(16-y)
plot(legend=false)
vectorfieldplot!((x,y) -> [1, F(y,x)], xlims=(x0, 10), ylims=(y0-1, y0+30))

g(x) =  16 - 13*exp(-0.2x)
plot!(g,  linecolor=:blue, linewidth=2)

correct me if I am wrong.