Julia ModelingToolkit: How to provide the initial conditions when an equation is provided in terms of a second derivative

I have the following code that works well and gives the correct result:

using ModelingToolkit
@parameters g vxᵢ
@variables t x(t) y(t) vy(t)
D = Differential(t)
eqs = [D(x) ~ vxᵢ,
       D(y) ~ vy,
       D(vy) ~ -g]

@named de = ODESystem(eqs, t, [x,y,vy], [g,vxᵢ])
ode_f = ODEFunction(de)

### Use in DifferentialEquations.jl

using OrdinaryDiffEq

α = π/2 / 90 * 45
vᵢ = 10.0
vxᵢ = vᵢ * cos(α)
vyᵢ = vᵢ * sin(α)

u₀ = [0.0, 0.0, vyᵢ]
tspan = (0.0, 3.0)
p = [9.81, vxᵢ]
prob = ODEProblem(ode_f,u₀,tspan,p)
sol = solve(prob, Tsit5(), abstol = 1e-4, reltol = 1e-4, saveat=0.01)

x = [u[1] for u in sol.u]
y = [u[2] for u in sol.u]
using Plots
p = plot(x,y)
gui(p)
println("Press any key")
read(stdin, Char)

I would like to change the equations as follows:

using ModelingToolkit
@parameters g vxᵢ vyᵢ
@variables t x(t) y(t)
D = Differential(t)
eqs = [D(x) ~ vxᵢ,
       D(vy)^2 ~ -g,
       # How to specify D(vy) at t=0 is vyᵢ
       ]

@named de = ODESystem(eqs, t, [x,y], [g,vxᵢ,vyᵢ])
ode_f = ODEFunction(de)

### Use in DifferentialEquations.jl

using OrdinaryDiffEq

α = π/2 / 90 * 45
vᵢ = 10.0
vxᵢ = vᵢ * cos(α)
vyᵢ = vᵢ * sin(α)

u₀ = [0.0, 0.0]
tspan = (0.0, 3.0)
p = [9.81, vxᵢ, vyᵢ]
prob = ODEProblem(ode_f,u₀,tspan,p)
sol = solve(prob, Tsit5(), abstol = 1e-4, reltol = 1e-4, saveat=0.01)

x = [u[1] for u in sol.u]
y = [u[2] for u in sol.u]
using Plots
p = plot(x,y)
gui(p)
println("Press any key")
read(stdin, Char)

How to specify D(y) at t = 0 is vyᵢ?

Don't directly create the ODEFunction. Follow the documentation and use the value maps. For example:

using ModelingToolkit, OrdinaryDiffEq
@parameters g vxᵢ vyᵢ
@variables t x(t) y(t)
D = Differential(t)
eqs = [D(x) ~ vxᵢ,
       (D^2)(y) ~ -g
]

@named de = ODESystem(eqs, t, [x, y], [g, vxᵢ, vyᵢ])
sde = structural_simplify(de)
tspan = (0.0, 3.0)

α = π / 2 / 90 * 45
vᵢ = 10.0
p = [g => 9.81, vxᵢ => vᵢ * cos(α), vyᵢ => vᵢ * sin(α)]
prob = ODEProblem(sde, [x => 0.0, y => 0.0, D(y) => vyᵢ], tspan, p)

sol = solve(prob, Tsit5(), abstol = 1e-4, reltol = 1e-4, saveat = 0.01)

solx = sol[x]
soly = sol[y]

# Better way to plot that
using Plots
plot(sol,vars=(x,y))
savefig("plot.png")