Exercise 10.1.5 with BVP

state as the original state and the first derivative

plt = plot(xaxis = L"x", yaxis = L"u(x)",
           title = "Different epsilons", leg=:bottomright);

domain = (eps(),1.0); # note that 1.0 vs 1 above for BVProblem vs shoot

for epsilon in [0.2, 0.02, 0.002]
  function ode!(f, u, epsilon, x)
    f[1] = -u[2]*(1-x) # first derivative
    f[2] = (u[1]^3 - u[1])/epsilon # second derivative
  end
  
  function bc!(g,u,epsilon,x)
    g[1] = u(0)[1] + 1.0
    g[2] = u(1)[1] - 1.0
  end
    
  bvp = BVProblem(ode!, bc!, [-1.0, 0.0], domain, epsilon);
  u = solve(bvp);
  plot!(u, idxs = 1, label="eps = $epsilon") # how to access it correctly so only first component is plotted? idxs = 1
end

plt