g = 9.81 ; % m/s/s force of gravity
m = 0.05 ; % mass of rocket
f = 16 ; % N upeard forece

%find velocity and height of the rocket in the first .15s, this loop starts when time in seconds is 0 and repeats the commands until 15 seconds 

t1 = 0;  % t1 = initial time
dt = 0.01; %dt = time increments
t = 0; v = 0; h = 0;  % to store values for plotting
i = 1; % counter
while t1<0.15  % First while loop. We loop for the first .15 secs
    a = f/m - g; % work out acceleration upwards
    v1 = a*t1;  % initial velocity
    h1 = .5*a*t1^2; % initial height
    t1 = t1 +dt;  % increment time by time increment
    i = i+1;  % increment counter
    t(i) = t1; v(i) = v1; h(i) = h1; % store time, vel and height in relevant arrays
    
end 

%this while loop models the velocity and height after 0.15 seconds until the parachute opens 

v2 = v1; % v2 = last height in stage one
t2 = t1; % time2 = last time from stage one
while v2 > -20  % Keep looping while vel > -20 Remember velocity is negative on the way down
    h2 = h1 + v1*(t2-t1)-.5*g*((t2-t1)^2); % Equation of motion s = ut +0.5atsquared
    v2 = v1 - g*(t2-t1);    % Equation of motion v = u + at
    t2 = t2 + dt;   % increment time by dt
    i = i+1;   % increment counter
    t(i) = t2; v(i) = v2; h(i) = h2; % store time, vel and height in arrays
end

%this while loop models the velocity and height from when the parachute opens till the rocket hits the ground. 
% Stage 3 of rocket (constant 20 m/s)

h3 = h2;  % h3 = last height in stage 2
t3 = t2;  % time3 = last time in stage 2
v3 = v2;  % v3 = last vel in stage 2
while h3>0  % keep looping whilst height > 0
    h3 = h2 + v3*(t3-t2); % equation of motion v = u + at
    t3 = t3+dt; % increment time
    i = i+1;  % increment counter
    t(i) = t3; v(i) = v3; h(i) = h3; % store time vel and height in array
end 

% Now to the plot

f = figure('visible','off');
subplot(2,1,1);  % The plot has two rows and one col. In placeholder one
plot(t, h),grid; % do this plot
xlabel('t(secs)'),ylabel('h(m)') % label the axes with units
title(' plot of height versus time') % plot y vs x
subplot(2,1,2) % in palceholder two do this plot
plot(t, v),grid; % plot vel vs time
xlabel('t'),ylabel('v') % label axes with units
title('plot of velocity versus time') 
saveas(f, 'rocket_plot.png');