clear
close all
T = [0 30]; % Solution range
g = 9.81; % Gravivt. Sign of g is built into model.
Y = [1000 100]; % Initial conditions [S V]
% Solve the system
[t,y] = ode45(@(t,y)sys1(t,y,g),T,Y);
%
% t is a vector of time values
% y is a matrix with 2 columns
% y(:,1) = column 1 = position
% y(:,2) = column 2 = velocity
% Plot results
figure(1)
hold on
plot(t,y(:,1),'b-')
plot(t,y(:,2),'r-')
hold off
legend('S','V')
xlabel('t')
% Obtain time when mass hits the ground by finding the
% time when the position is zero
yval = 0;
tground = revinterp(t,y(:,1),yval)
% Obtain the time when the mass reaches is maximum height
% by finding the time when the velocity is zero.
vval = 0;
[tmax,tind] = revinterp(t,y(:,2),vval)
% Obtain the maximum height using quadratic interpolation
% into the position array. Substitute tmax into the
% resulting quadratic function.
%
% The last 2 in the calling sequence indicates that we want
% a second order (i.e., quadratic) polynomial fit.
p = polyfit(t(tind:tind+2),y(tind:tind+2,1), 2);
smax = polyval(p,tmax)
% Maximum height obtained by doing simple linear interpolation
% into the position array. Notice that this is less accurrate
% than the quadratic value.
smax2 = revinterp(y(:,1),t,tmax)