% Class 26
% Example 24.21 after class
%
% Step 1) manipulate equation so it is two first order ODE
% dx/dt = h, and dh/dt = -c/m*h+g
%
% Step 2) make sub function for ODE solver (see SYS below)
% Step 3) make sub function to find inital value of slope (slopefind)
% Step 4) use FZERO to solve for the unknown slope
slope = fzero(@slopefind,1);
% Step 5) use the slope to solve the probelm (e.g. plot how x varries with
% time!)
[tp,yp] = ode45(@SYS,[0 12],[0 slope]);
plot(tp,yp(:,1))
xlabel('Time (s)')
ylabel('Distance (m)')
%
% NOTE: This equation is linear, so two shots with interpolation 
% will also work.  Check that idea here!
% First guess
slope1 = 10;
[~,yp] = ode45(@SYS,[0 12],[0 slope1]);
Ans1 = yp(end,1);
% Second guess
slope2 = 20;
[~,yp] = ode45(@SYS,[0 12],[0 slope2]);
Ans2 = yp(end,1);
% Interpolate to get real slope
x = [Ans1 Ans2]'; y = [slope1 slope2]';
slope_real = interp1(x,y,500,'linear','extrap'); % NOTE, need to include this for interp1 to extrapolate beyond bounds
% Now solve the problem!
[tp,yp] = ode45(@SYS,[0 12],[0 slope_real]);
hold on
plot(tp,yp(:,1),'+')

function E = slopefind(sg)
[tp,yp] = ode45(@SYS,[0 12],[0 sg]);
result = yp(end,1);
desired = 500;
E = result - desired;
end

function dy = SYS(t,y)
% t = t;
% y(1) = x
% y(2) = h
c = 12.5; %kg/s
m = 70; %kg
g = 9.91; %m/s^2
dy = [y(2);...
    -c/m*y(2)+g];
end