% Class 28
% Coded by NFR on 4.25.19
close all
clc
% Warm up
[tp,yp] = ode23s(@SYS,[0 6000],[1 1]);
%xlabel('Time')
ylabel('y_1')
%
% Question 24.3 from book
% Example of LINEAR boundary value problem
% Solve first guess point
sg1 = 0.1;
[tp,yp] = ode45(@SYS2,[0 20],[5 sg1]);
gp1 = yp(end,1)
% Solve second guess point
sg2 = 1;
[tp,yp] = ode45(@SYS2,[0 20],[5 sg2]);
gp2 = yp(end,1)
% Use linear interpolation of these two points
% to determine the unknown slope
Ratio = (sg1-sg2)/(gp1-gp2);
Slope = sg1-Ratio*(gp1-8);
% Solve the problem with the interpolated slope
[tp,yp] = ode45(@SYS2,[0 20],[5 Slope]);
result = yp(end,1)
% Pause on the interpolation method
% Solving as root problem
Slope = fminsearch(@OPT,1)
% Test this slope
[tp,yp] = ode45(@SYS2,[0 20],[5 Slope]);
result = yp(end,1)
plot(tp,yp(:,1))
xlabel('x')
ylabel('y')
% Second example of BVP
% 24.16 from the textbook
% First step - create a subfunction for the ODE
% Second step - create a subfunction to solve for unknown slope
% Third step - use unknown slope to solve ODE
% Fourth - check the end value
% Fifth - plot the solution
slope = fzero(@OPT2,5);
[tp,yp] = ode45(@SYS3,[0 3],[0 slope]);
plot(tp,yp(:,1))

function E = OPT2(sg)
[tp,yp] = ode45(@SYS3,[0 3],[0 sg]);
result = yp(end,1);
desired = 0;
E = (result-desired);
end

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



function E = OPT(sg)
[tp,yp] = ode45(@SYS2,[0 20],[5 sg]);
result = yp(end,1);
desired = 8;
E = (result-desired)^2;
end

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


function dy = SYS(t,y)
% t = t
% y(1) = y
% y(2) = h
mu = 1000;
dy = [y(2);...
    mu*(1-y(1)^2)*y(2)-y(1)];
end