% Uppgift 9

%% Shooting method
% The "shooting method" is a scattershot approach to solving ODEs that are difficult to solve more analytically. It reduces a boundary value problem (solution that satisfies boundaries) to an initial value problem (which initial value satisifies the boundary?).
% Practically, we solve for many initial values s in a "scattershot" manner and then build up a set of endpoints q for each initial value.
% By finding the roots of q, one finds which initial values satisfy the endpoint boundary.
% Lös randvärdesproblemet
% {u^{\prime } }^{\prime } +u^2 =1
% u\left(0\right)=u\left(1\right)=0
% Skriv om som system av ODE
% {u_1 }^{\prime } =u_2
% {u_2 }^{\prime } =1-u_1^2
% u_1 \left(0\right)=0,\;u_2 =s
%
% This problem is non-linear and non-homogenous, making it tricky to solve using the finite-differences method.
% OBS! This file is messy because matlab makes it hard to pass functions around.

close all

s = linspace(-20,50)
f1 = @(u1,u2) u2;
f2 = @(u1, u2) 1-u1.^2;
f = @(t,U) [f1(U(1),U(2)); f2(U(1),U(2))];

tspan = [0 1]; % The bounds.
q = zeros(size(s));
for k = 1:length(s)
    [t u] = ode45(f, tspan, [0; s(k)]);
    q(k) = u(end,1)-0;
end

plot(s, q)
xlabel('s')
ylabel('q(s)')
grid on
ylim([-10 10])

% Clearly, two solutions for q(s) satisfy the boundary conditions.
q_root1 = fzero(@(s) findq(s), [-1 0] )
q_root2 = fzero(@(s) findq(s), [32 34] ) 
hold on
plot(q_root1, 0, 'ro', q_root2, 0, 'ro' )


[t u] = ode45(@odesystem_uppg9, tspan, [0 q_root1]);
figure
plot(t,u(:,1))
hold on
%plot(t, u(:,2))
grid on
title(['s = ' num2str(q_root1)])


[t u] = ode45(@odesystem_uppg9, tspan, [0 q_root2]);
plot(t,u(:,1))
hold on
%plot(t, u(:,2))
grid on
title(['s_1 = ' num2str(q_root1) ' s_2 = ' num2str(q_root2)])
xlabel('t')
ylabel('u(t)')

%
%% Finite differences method
% This differential equation is non-linear and so the usual method of Au=b won't suffice.
% TODO: solve with relaxation method http://www.ohiouniversityfaculty.com/youngt/IntNumMeth/lecture34.pdf
%

g0 = 0;
g1 = 0;
L = 1;
N = 30;

%

function q = findq(s)
[t u] = ode45(@odesystem_uppg9, [0 1], [0; s]);
q = u(end,1) - 0;
end

function uprim = odesystem_uppg9(t,U) 
    uprim = [U(2); 1-U(1).^2];
end