% Second Order Ordinary Differential Equations - Initial Value Problems
% Consider the equation
% \frac{d^2 \mathrm{y}}{\mathrm{d}t^2 }=-g+\frac{Cv^2 }{m}
% with the inital conditions y(0)=h and \dot{y} \left(0\right)=0. The built-in function ode45 will also solve this equation, but first we have to put it in the form of a system of first order differential equations. The simplest way to do so is to introduce an intermediate variable, z, defined as
% \frac{\textrm{dy}}{\textrm{dt}}=z
% Then we can write
% \frac{d^2 y}{dt^2 }=\frac{\textrm{dz}}{\textrm{dt}}=-g+\frac{Cv^2 }{m}
% and our system of equations becomes
% \begin{array}{l}
% \frac{\textrm{dy}}{\textrm{dt}}=z\\
% \frac{\textrm{dz}}{\textrm{dt}}=-g+\frac{Cv^2 }{m}
% \end{array}
% with initial conditions
% \begin{array}{l}
% y\left(0\right)=h\\
% z\left(0\right)=0
% \end{array}
% Now we can solve this system as follows

h=600;
tr=[0 3];
initv=[h 0];
[t,y]=ode45(@rkfalling, tr, initv);
plot(t,y(:,1), 'LineWidth', 2)
set(gca,'FontSize',20)
ylabel('y')
xlabel('time')
figure
plot(t,y(:,2), 'LineWidth', 2)
set(gca,'FontSize',20)
ylabel('dy/dt')
xlabel('time')

%
% The trange vector indicates that we will solve the equation on the range 0<t<3. The init vector sets the initial values of both y and z.
%
% As called, ode45 expects to find the function definition in a user-defined function which I've called rkfalling. For this problem, rkfalling is as follows:

function r=rkfalling(t,w)
mass=80;
C=20;
g=9.81;
r=zeros(2,1);
r(1)=w(2);
r(2)= C*w(2).^2/mass-g;
end

% Here w is a vector holding a value for both y and z. I've chosen y to be w(1) and z to be y(2). If this is the order of the variables in w, then rkfalling must return a vector which I've called r which holds values for dy/dt and dz/dt (in that order).
% © James Blanchard, University of Wisconsin-Madison, blanchard@engr.wisc.edu