Although we can perform many modeling tasks in Excel, a full$-$fledged scientific computing

programming language like MATLAB can do much more and do so more efficiently. In this

problem, we will learn how to write simple MATLAB programs to solve an ODE numerically.

$($a$)$ In the last module, we learned how to use Euler method to solve an ODE numerically:

$\frac{dy}{dt}=f(t,y),$;$,y(t=0)=y_0$

Write a MATLAB function that implements Euler method. The definition should be:

$[tt,yy]=$ feuler$(fun,y0,tf,h)$

where fun is a user$-$defined function that should take in two arguments, $t$ and $y,$ in that

order, and return $f(t,y).$ The program will simulate the ODE from $t=0$ to $t=t_f$ with

time step $h,$ and return the simulated data points $(t_i,y_i)$ as two vectors tt and yY$,$

consisting of the $t_i\text{'}$s and the corresponding $y_i\text{'}$s$,$ respectively.

$($b$)$ The Euler method we learned is properly called the "forward Euler" method since it

assumes that the slope between $(t_i,y_i)$ and $(t_{i+1},y_{i+1})$ can be approximated by the slope

evaluated at $(t_i,y_i).$ The formula is:

$y_{i+1}=y_i+hf(t_i,y_i)$

Another variant of Euler method is the "backward Euler" method, which uses the slope

evaluated at the point $(t_{i+1},y_{i+1})$ instead, resulting in the formula:

$y_{i+1}=y_i+hf(t_{i+1},y_{i+1})$

Since we actually do not know $y_{i+1}$ yet, we cannot evaluate $f(t_{i+1},y_{i+1})$ right away and

plug it into the formula. Instead, we need to solve this implicit formula to calculate $y_{i+1}.$

In MATLAB, you can do this with the fzero function $($see the MATLAB tutorial video on

"Root Finding"$).$

Write a MATLAB function that implements "backward Euler" method. The definition

should be:

$[$tt$,$ yy$]=$ beuler$($fun$,$ y$0,$ tf$,$ h$)$

$($c$)$ Yet another method for solving ODEs is the "mid$-$point method," which uses the slope at

the mid$-$point between $t_i$ and $t_{i+1}$ :

$y_{i+1}=y_i+hf(\frac{t_i+t_{i+1}}{2},\frac{y_i+y_{i+1}}{2})$

Similar to the backward Euler method, you will need to use fzero to find $y_{i+1}$ at each

step. Write a MATLAB function that implements this method. The definition should be:

$[$tt$,$ yy$]=$ meuler$($fun$,$ y$0,$ tf$,$ h$)$

$($d$)$ Apply each of the methods from Parts $($a$)-($c$)$ to Question $1($c$).$ Write a MATLAB program

to solve the ODE of Question $1($c$)$ for $C(t)$ using each method one by one, and plot the

solutions from each method, together with the analytical solution in the same plot. For

simplicity, calculate $V(t)$ directly from its analytical solution and only use the numerical

methods to solve for $C(t).$ Your function definition should be:

$[]=tank(Cin,$ Fin, $Vo,k,tf,h)$

$($e$)$ For a parameter setting of your choice, use your program to compare and discuss the

error behaviors of the $3$ methods. Show several plots of the solutions at different step

sizes $h.$ Make sure you try a large enough $h$ so that you can see the differences more

clearly.