Implement three Julia functions for numerical solution of 1st-order IVPs

\begin{align*} y' &= f(t,y) \\ y(t_0) & = y_0 \end{align*}for forward Euler, Runge-Kutta 2nd order (midpoint method), and Runge-Kutta 4th order. The functions should take inputs `f, y0, t0, Δt, N`

and returns $N+1$-dimensional vectors `t`

and `y`

, where `t`

is the vector of discrete times $t_n = t_0 + n \Delta t$ and `y`

is the vector of approximate values of $y$ at those times,
$y(t_0), y(t_1), y(t_2), \ldots$.

For `forwardeuler`

and `rungekutta2`

, feel free to copy code from the lecture notebook.

In [ ]:

```
function forwardeuler(f, y0, t0, Δt, N)
end
function rungekutta2(f, y0, t0, Δt, N) # a.k.a. midpoint method
end
function rungekutta4(f, y0, t0, Δt, N)
end
```

Compare the three numerical solution methods on the initial value problem

\begin{align*} y' &= t y + t^3 \\ y(0) &= 1 \end{align*}to the true solution $y(t) = [2 + y(0)] \, e^{t^2/2} - t^2 - 2$ over the interval $0 \leq t \leq 2$. Your comparison should consist of two plots. The first plot should show $y(t)$ versus $t$ for the true solution and the three numerical solutions. The second plot should show the global errors of the three methods at the last time step versus the time step $\Delta t$.

For guidance on the plots and the structure of the computations to produce them, look at the lecture notebook for 2016-12-02.

In [ ]:

```
```

Revise your three numerical IVP solvers to work when $y$ and $f$ are vector-valued. If $d$ is the dimension of $y$ and $f$, the `y`

output should now be an $d \times (N+1)$ matrix, whose columns are the $d$-dimensional vectors $y(t_0), y(t_1), y(t_2), \ldots$.

In [ ]:

```
```

The equations of motion for a small satellite orbiting a much more massive central body are

\begin{align*} x'' &= -\frac{Gmx}{(x^2+y^2)^{3/2}} \\ y'' &= -\frac{Gmy}{(x^2+y^2)^{3/2}} \end{align*}In these equations, $(x,y)$ is the position of the satellite in the plane of the orbit, $G$ is the universal gravitational constant, and $m$ is the mass of the central body at the origin.

**(a)** Convert this system of two second-order ODEs to a system of four first-order ODEs in the four variables $x, y, v_x, v_y$, where $v_x = x'$ and $v_y = y'$.

**(b)** Write a Julia function `f(t,x)`

that maps the vector $x = [x, y, v_x, v_y]$ into $x' = [x', y', v_x', v_y']$
according to the system of four first-order ODEs you derived in **(a)**. Since this is a math class and not a physics class, you can set the physical constant $G$ and $m$ to 1.

**(c)** Compute numerical solutions to the system of ODEs using the initial condition

using Forward Euler, 2nd-order Runge-Kutta, and 4th-order Runge-Kutta, using $\Delta t = 0.04$ over the interval $0 \leq t \leq 10$

**(d)** Make a plot of the $x(t), y(t)$ orbit of the satellite as computed using the three time-stepping algorithms. Plot Forward Euler with a blue line, 2nd-order Runge-Kutta in red, and 4th-order Runge-Kutta in green. Make sure the axes are equispaced in $x$ and $y$ by using `axis("equal")`

, and crop the plot to $-1 \leq x \leq 1.5$ and $-1.5 \leq y \leq 1$ using `xlim(-1, 1.5); ylim(-1.5, 1)`

.

**(e)** Comment on your results. Describe the differences between the three approximate orbits. Which of them makes most sense from a physical perspective, and why?

In [ ]:

```
```