% Class 26 
% ODE Boundary Value Problems
%
% Shooting method with the BUICK
x = [23.6 11.4];
y = [67.8 80.1];
% Determine the desired angle using interpolation
yd = interp1(x,y,17.49)
%
% Use shooting method to solve 24.3
% 1) Converted second order equation into 
% two first order
% 2) Created a system of equations (local function)
% 3) Solve a first guess point
hg1 = 10;
[tp,yp] = ode45(@SYS1,[0 20],[5 hg1])
ans1 = yp(end,1);
% 4) To interpoalte solve a second guess point
hg2 = 15;
[tp,yp] = ode45(@SYS1,[0 20],[5 hg2])
ans2 = yp(end,1);
% 5) Interpolate using these two solution to 
% find the actual value of h that allows the 
% solution to converge to end boundary cond.
x = [ans1 ans2]
y = [hg1 hg2]
% Querry point should be the value at the end boundary
hact = interp1(x,y,8,'linear','extrap')
% Now we know what initial condtion to use 
% for h to get the solution, solve the problem!
[tp,yp] = ode45(@SYS1,[0 20],[5 hact]);
% PLOT to see if it makes sense
hold off
plot(tp,yp(:,1))
xlabel('x')
ylabel('y')
% PAUSE on the interpolation
% Solve this problem using root finding
% We need to set up another local funtion
% that we can call in a root finding solution
hact = fzero(@OPT1,1)
% Use this to solve the problem
[tp,yp] = ode45(@SYS1,[0 20],[5 hact]);
plot(tp,yp(:,1))
xlabel('x')
ylabel('y')
%
% Example 24.21 - code will be posted (similar to above)
%
% Example 24.16
% 1) Transformed second order into two first order
% 2) Make system of equations
% 3) set up our root finding local function
% 4) guess the slope (h) value, and use root fiding
hg = 1;
hact = fzero(@OPT2,hg);
% 5) Use this actual slope value at the start
% to now solve the problem
[tp,yp] = ode45(@SYS2,[0 3],[0 hact]);
% 6) Visualize the solution with plot
plot(tp,yp(:,1)) % Again, your y vector is the first column
% the second column are the h values (slopes) that we don't care about
% 7) If you have an analytical solution
% you can also plot that on top to see
% how good your numerical solution is!

function E = OPT2(hg)
[tp,yp] = ode45(@SYS2,[0 3],[0 hg])
endvalue = yp(end,1);
desiredvalue = 0;
E = endvalue - desiredvalue;
end

function dydt = SYS2(t,y)
% y(1) = y
% y(2) = h
% t = x
E = 200*10^9; %Pa
I = 30000/(100^4); %m^4
w = 15000; %N/m
L = 3; %m
dydt = [y(2);
    (w*L*t/2-w*t^2/2)/(E*I)];
end

function E = OPT1(hg)
[tp,yp] = ode45(@SYS1,[0 20],[5 hg]);
EndPoint = yp(end,1);
Desired = 8;
E = EndPoint-Desired;
end

function dydt = SYS1(t,y)
% y(1) = y
% y(2) = h
% t = x
dydt = [y(2);
    (2*y(2)+y(1)-t)/7];
end