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)