%% LAB 2
% Table of Contents
% 2. Initial-value Problems for 1st-order Differential Equations
%     2.1 Riktningsfält
%     2.2 Differensmetoder
%     2.3 Konvergens
%     Uppgift 2
%     Uppgift 3: Jämför euler framåtmetod med ode45().
% 2.5 System av ODE med ode45()
%     Fasporträtt: u_1 vs. u_2
%     Limit cycle
%     Uppgift 4: same procedure for new system of ODE's.
% 2.6: Högre ordningens diferentialekvationer
%     Uppgift 6: Damped pendulum
% Randvärdesproblem/Boundary-value problems
%     3.2 Inskjutningsmetoden/Shooting method (https://en.wikipedia.org/wiki/Shooting_method)
%         Solve the previous heat transfer problem using the Shooting method
%% 2. Initial-value Problems for 1st-order Differential Equations
% Initial-value problems for first order differential equations have the form
% u^{\prime } =f\left(t,u\right)
% u\left(0\right)=u_0, where u0 is a constant.
% What makes it a differential equation is that u' is a function of u. The function will evolve differently depending on the initial value.
%% 2.1 Riktningsfält
% As an example, we consider the 1st-order DE u^{\prime } \left(t\right)=-u\left(t\right)+\sin \left(t\right)+\cos \left(t\right), u\left(0\right)=u_0.
% We have a derivative defined in terms of

f = @(t,u) -u+sin(t)+cos(t); % högerled. dvs. u'(t)
f_exakt=@(t,u0) sin(t)+u0.*exp(-t); % Analytic uprim.

% Draw vector-field of \frac{\textrm{du}}{\textrm{dt}}. This shows how the function will change at each point.
% That is, the axes are u(t) and t.

close all

a=0;
b=4;
N=20;
t=linspace(a,b,N); 
u=linspace(-2,2,N); % u(t)
[T,U]=meshgrid(t,u); % Grid of (t, u(t))
DT=ones(size(T)); % d/dt(t) = 1
DU=f(T,U); % For our range of t and u(t), find u'(t).
figure
quiver(T,U,DT,DU,0.9) % Plot vector [dt du] at each point.
hold on
xlabel('t');
ylabel('u(t)')
hold on
grid on

plot(t, f_exakt(t, 0))
plot(t, f_exakt(t, 1))
plot(t, f_exakt(t, 2))
plot(t, f_exakt(t, -2))

%
%% 2.2 Differensmetoder
% We approximate the derivative at a point using a forward difference quotient: a linearisation between two adjacent elements.
% D_+ u\left(t_n \right)=\frac{u\left(t_{n+1} \right)-u\left(t_n \right)}{h}\approx u^{\prime } \left(t_n \right)=f\left(t_n ,u\left(t_n \right)\right), an approximation of the definition given at the top of the script.
% Rearranging gives u\left(t_{n+1} \right)\approx u\left(t_n \right)+h\;f\left(t_n ,u\left(t_n \right)\right)
% If we let u_nbe an approximation of u\left(t_n \right), then Euler's forward-method is
% u_{n+1} =u_n +h\;f\left(t_n ,u_n \right)
% This method is explicit since all variables on the r.h.s. are known.
%
% We derive using Euler's forward-method, using initial value u_0, and plot the curve.

close all

f =       @(t,u) -u+sin(t)+cos(t); % u'(t) = f(t, u)
f_exakt = @(t,u0) sin(t)+u0.*exp(-t);

N = 15;
h = (b-a)/N;
t = linspace(a,b,N+1);
u = zeros(size(t));
u0 = 1;
u(1) = u0; % Initial value.

for n = 1:N
    u(n+1) = u(n)+h.*f(t(n), u(n));
end
plot(t, u);
hold on
xlabel('t');
ylabel('u(t)')
grid on

% Euler's backwards-method is derived in a similar way:
% u_{n+1} =u_n +h\;f\left(t_{n+1} ,u_{n+1} \right)
% This method is implicit because u_{n+1}is unknown and appears on both sides!
% In each step, a non-linear equation g\left(w\right)=w-u_n -h\;f\left(t_{n+1} ,w\right)must be solved using eg. newton's method.
%
% The trapezium method combines Euler's backwards and forwards methods, and is thus also implicit.
% u_{n+1} =u_n +\frac{h}{2}\left(f\left(t_n ,u_n \right)+f\left(t_{n+1} ,u_{n+1} \right)\right)
%
% Heun's Method takes the trapezium method and substitutes the implicit u_{n+1}for a forward method approximation (and substitutes t_{n+1} \to t_n +h), making it explicit.
% u_{n+1} =u_n +\frac{h}{2}\left(f\left(t_n ,u_n \right)\;+f\left(t_n +h,u_n +h\;f\left(t_n ,u_n \right)\right)\right)
%
% Finally, a brute force method is to simply integrate the differential equation for each time step:
% u\left(t_{n+1} \right)=u\left(t_n \right)+\int_{t_n }^{t_{n+1} } u^{\prime } \left(t,u\left(t_n \right)\right)\textrm{dt}

u = zeros(size(t));
u(1) = u0;
for n = 1:N
    u(n+1) = u(n) + h/2*( f(t(n),u(n)) + f(t(n)+h, u(n)+h*f(t(n),u(n))) ); % Heun
end
plot(t, u)

title(['u_0 = ', num2str(u0)], ' ')

integral_method = zeros(size(t));
integral_method(1) = u0;
for n = 1:N
    integral_method(n+1) = u(n) + integral(@(t) f(t,u(n)), t(n), t(n+1)); % Integrate the differential equation in each each time-step.
end
plot(t,integral_method)

plot(t,f_exakt(t,u0))

%plot(t, f(t, 0), 'g--',t, f(t, -1), 'g--',t, f(t, 1), 'r-',t, f(t, 2), 'g--'  )
%hold off
hold off
legend("Euler forward", "Heun", 'Integration using integral()', 'Exact')

%% 2.3 Konvergens
% All the difference methods described above converge, and thus we can get arbitrarily close to the exact value as h\to 0.
% Euler's methods have linear convergence, while Heun's method has quadratic convergence.
% That is, if we halve the step-size, the error will be 1/2 and 1/4, respectively.
%
%% Uppgift 2
% Lös problemet med Eulers framåtekvation. Välj en bra steglängd. Rita riktningsfältet.

close all
uprim = @(t,u) cos(3.*t)-sin(5*t).*u;

a = 0;
b = 5;
N = 30;
t = linspace(a,b,N); % För N intervall krävs N+1 punkter.
u = linspace(0, 3, N);
[T,U]=meshgrid(t,u);
DT = ones(size(T)); 
DU = uprim(T,U);
quiver(T,U,DT,DU, 0.8), 
hold on

% Vi ser att Euler metoden konvergerar mot den exakta lösningen i takt med att steglängden minskas.
% Vi kan även se att felet ackumuleras med tid (dvs. med varje nytt steg i diff.ekvationen).

for N=[10 20 30 50 100 1000 10000]
    h=(b-a)/N;
    t=linspace(a,b,N+1);u=zeros(size(t));
    ua=2;
    u(1)=ua; % at u(t), t=0
    for n=1:N
        u(n+1)=u(n)+h*uprim(t(n), u(n));
    end
    plot(t, u)
end
grid on
xlabel('t');
ylabel('u(t)')
hold off
title('Effects of time-step-size h using Euler-Forward')

%
% We can also look at the resulting value of the first equation u=-u^{\prime } +\sin \;t+\cos \;t at t=4.
% Depending on the method and step-size used, different errors will be introduced.
% We will look at the difference between forward-Euler, Heun and the analytic solution at t=4.

close all

u0 = 1;
u_exakt = @(t,u0) sin(t) + u0*exp(-t);
uprim = @(t,u) -u +sin(t) + cos(t);
a = 0;
b = 4;
h_vec = [0.4 0.2 0.1 0.05 0.025 0.01 0.001];
N = (b-a)./h_vec;

res_eulerForward = zeros(size(N));
res_heun = zeros(size(N));
res_exakt = 0

for i = 1:length(N)
    t = linspace(a,b,N(i)+1);
    u = zeros(size(t));
    h = (b-a)./N(i);
    u(1) = u0;
    % Euler forward
    for n=1:N(i)
        u(n+1)=u(n)+h*uprim(t(n), u(n));
    end
    %plot(t, u); hold on
    res_eulerForward(i) = u(end);

    u = zeros(size(t));
    u(1) = u0;
    for n = 1:N(i)
        u(n+1) = u(n) + h/2*( uprim(t(n),u(n)) + uprim(t(n)+h, u(n)+h*uprim(t(n),u(n))) );
    end
    res_heun(i) = u(end);
end
error_eulerForward = res_eulerForward - u_exakt(4,u0);
plot(h_vec, error_eulerForward); hold on
error_heun = res_heun - u_exakt(4,u0);
plot(h_vec, error_heun)
grid on
title("Convergence of forward Euler and Heun")
xlabel("h")
set ( gca, 'xdir', 'reverse' ) % Flip x axis in plot.
legend('Euler Forward', 'Heun')

%
%% Uppgift 3: Jämför euler framåtmetod med ode45().

close all
figure
uprim=@(t,u) cos(3.*t)-sin(5*t).*u;
a=0;
b=5;
N=20;
t=linspace(a,b,N+1); % För N intervall krävs N+1 punkter.
u=linspace(0, 3, N);
[T,U]=meshgrid(t,u);
DT=ones(size(T)); DU=uprim(T,U);
quiver(T,U,DT,DU, 0.8), 
hold on

h=(b-a)/N;
u=zeros(size(t));
u(1)=ua; % at u(t), t=0
for n=1:N
    u(n+1)=u(n)+h*uprim(t(n), u(n));
end
plot(t, u), hold on

[t,u]=ode45(uprim,[a b], ua);
plot(t,u), grid on

legend('','Eulers framåtmetod', 'ode45')

%
%% 2.5 System av ODE med ode45()
% Previously, our definition of an ODE was u^{\prime } \left(t\right)=f\left(t,u\left(t\right)\right), a\le t\le b, u\left(a\right)=u_a.
% This can be extended to systems of ODE's without changing the definition, by using vector notation.
% \overrightarrow{u^{\prime } \left(t\right)} =\overrightarrow{f} \left(t,\overrightarrow{u} \left(t\right)\right), \overrightarrow{u} \left(a\right)=\overrightarrow{u_a }
% The above says:
%  \overrightarrow{u}is a system of equations \left\lbrack u_1 ,u_2 ,{\ldotp \ldotp \ldotp u}_m \right\rbrack that we are trying to obtain/numerically compute.
% \overrightarrow{f} \left(t,\overrightarrow{u} \right)is a system of differential equations that are functions of t and the equations in \overrightarrow{\;u}.
% \overrightarrow{u_a }are the initial values of the equations defined in \overrightarrow{u}.
%
% Ex.
% We have the following system of diff. equations:
% {u_1 }^{\prime } =2u_1 -u_1 u_2, u_1 \left(0\right)=0\ldotp 1
% {u_2 }^{\prime } =u_2 +0\ldotp 4u_1 u_2 -u_2^2, u_2 \left(0\right)=0\ldotp 1
% After using ode45(), we obtain numerical approximations of \overrightarrow{u}, the original functions.
% \overrightarrow{u} =\textrm{ode45}\left(\overrightarrow{f} \left(t,\overrightarrow{u} \right),t,\overrightarrow{u_0 } \right)

close all
f1 = @(u1,u2) 2*u1-u1.*u2; % The differential equations.
f2 = @(u1,u2) u2+0.4*u1.*u2-u2.^2;
f = @(t,u) [f1(u(1),u(2)); f2(u(1),u(2))]; % Vectorized. OBS! odefun must accept arguments (t,y), even if one of them is not used!
u0 = [0.1; 0.1]; % initial values.
T = 8; % Max time bound.
tspan = [0 T];
[t, U] = ode45(f, tspan, u0);
plot(t,U(:,1), 'b', t, U(:,2), 'g'), hold on, grid on
hold off
legend('u_1(t)', 'u_2(t)')

% This system of ODEs is autonomous, that is u^{\prime } \left(t\right)=f\left(u\left(t\right)\right), since the derivative does not directly depend on the indpendent variable t, only indirectly through u(t). That means that the evolution of the system depends only on the current state of the system.
% This means that stationary points can be found. The stationary points are shown in the phase portrait.
% If the independent variable is time, the system is said to be time-invariant.
%% Fasporträtt: u_1 vs. u_2
% The phase portrait can be interpreted as visualizing the trajectories of solutions to the system depending on initial conditions for \overrightarrow{u}.
% In the phase portrait below, we see three stationary points. One is stable, corresponding to the final, stable solution of the time-plot above. Some initial conditions will evolve to this stable point, others will diverge.

u1 = linspace(-1,4,25); 
u2 = linspace(-1,3,25);
[U1,U2] = meshgrid(u1,u2);
quiver(U1,U2,f1(U1,U2),f2(U1,U2), 1.5), hold on, grid on % vector components are the derivatives u1' and u2'.
plot(U(:,1), U(:,2), 'r') % Plot the solution that corresponds to a particular initial condition u0.

% To illustrate, change initial condition and recompute the system of ODEs.
% A different trajectory will be followed. 
% It stabilizes at the same point (sink).
 u0 = [0.1 0.5 -0.5 2; 0.5 3 3 3]; % Iterate over these pairs of points (columns are the pairs).
for i = 1:length(u0)
[t, U] = ode45(f, tspan, u0(:,i));
plot(U(:,1), U(:,2), 'g-', 'LineWidth',2)
end
xlabel('u_1'); ylabel('u_2');
axis([-1 4 -1 3])
plot(0,0,'r.', 0,1,'r.', 2.5 ,2,'r.')
legend('', 'u_0=[0.1;0.1]', 'u_0=[0.1;0.5]')

%
%% Limit cycle
% As another example, we solve another autonomous that results in a single unstable stationary point.
% Solutions (based on initial conditions) will follow a trajectory that spirals into a limit cycle (a closed trajectory that other solutions spiral into as t\to \infty).
% The stationary point resides in the interior of the limit cycle.

close all
f1 = @(u1, u2) -u1 + 0.04*u2 + u1.^2.*u2;
f2 = @(u1, u2) 0.8 - 0.04*u2 - u1.^2.*u2;
f = @(t,u) [f1(u(1),u(2)); f2(u(1),u(2))]; % System of diff.equations.

u1 = linspace(0,4.5,20);
u2 = linspace(0, 3, 20);
[U1, U2] = meshgrid(u1, u2);
quiver(U1, U2, f1(U1, U2), f2(U1, U2),5.5) % How to normalize arrow length? ie. scale each arrow 
hold on

tspan = [0 12];
u0 = [0.1 0.5 1 2 1 ; 0.5 3 3 3 1.5]; % Iterate over these pairs of points (columns are the pairs).
for i = 1:length(u0)
    [t, U] = ode45(f, tspan, u0(:,i));
    plot(u0(1,i), u0(2,i), 'ko')
    plot(U(:,1), U(:,2), 'LineWidth',2)
end
grid on
plot(0.8,1.1765, 'r.', 'LineWidth', 2)
axis([0 4.5 0 3])

%
%% Uppgift 4: same procedure for new system of ODE's.
% This system has two stable stationary points.

close all
% Funktioner u'1 och u'2 i ekvationssystem
f1=@(u1,u2) cos(1+u1-u2.^2)+0.3;
f2=@(u1,u2) sin(0.3+u1+u2)-0.4*u1;
f=@(t,u)    [f1(u(1),u(2)); f2(u(1),u(2))];
% Begynnelsevärden u1(0) och u2(0)
u0=[1.5; 3];
% Lösningsintervall
T=8; tspan=linspace(0,T,200);
% Lös ekv.systemet med ode45 på interval tspan med begynnelsevärden u0.
[t,U]=ode45(f,tspan,u0);
% Tidsfunktion
subplot(1,3,1)
plot(t,U(:,1),t,U(:,2)), grid on, hold on, legend('u_1', 'u_2'), xlabel('t')
% Riktiningsfält
% First, set up solution space
u1=linspace(0,3,30);
u2=linspace(0,3,30);
[U1,U2]=meshgrid(u1,u2);
% Then use quiver with the individual functions, and our grid as points, as
% parameters.
subplot(1,3,2:3)
quiver(U1,U2,f1(U1,U2),f2(U1,U2), 0.9), hold on
axis([0 3 0 3]), axis square, xlabel('u_1'), ylabel('u_2');

% Riktningsbanor
for	n1 = 0:1:3
    for n2 = 0:1:3
        u0(1)=n1;
        u0(2)=n2;
        [t,U]=ode45(f,tspan,u0);
        plot(U(:,1),U(:,2), 'r', 'LineWidth', 1.5)
    end
end
grid on

%
%% 2.6: Högre ordningens diferentialekvationer
%
% We can rewrite a higher order diff.equation as a system of first-order diff.equations.
% If we let:
% u_1 =u,
% u_2 =u^{\prime },
% u_3 ={u^{\prime } }^{\prime }
% and write the above a system of ODEs:
% \overrightarrow{u^{\prime } } =\overrightarrow{f} \left(t,\overrightarrow{u} \right)
% \overrightarrow{u} =\left\lbrack \begin{array}{c}
% u_1 \\
% \ldotp \ldotp \ldotp \\
% u_n
% \end{array}\right\rbrack, \overrightarrow{f} \left(t,\overrightarrow{u} \right)=\left\lbrack \begin{array}{c}
% u_2 \\
% \ldotp \ldotp \ldotp \\
% f\left(t,u_1 ,\;\ldotp \ldotp \ldotp ,u_n \right)
% \end{array}\right\rbrack
% The system will have the same number of equations as the order to the diff.equation.
%
% We demonstrate by solving the "mathematical" pendulum (ie. ideal -> no losses -> perpetual motion!)
% For a pendulum of mass m, length L from pivot, at angle \phi, under influence of gravity g.
% \textrm{mL}\ddot{\varphi} \left(t\right)=-\textrm{mg}\;\sin \left(\varphi \left(t\right)\right)
% Starting angle \varphi \left(0\right)=\varphi_0
% Starting angular velocity \omega \left(0\right)=\dot{\varphi} \left(0\right)=0
%
% Let u_1 =\varphi ,u_2 =\varphi^{\prime } ,u_3 ={\varphi^{\prime } }^{\prime }
% {u_1 }^{\prime } =u_2 ,{\;\;\;\;u}_1 \left(0\right)=\varphi_0
% {u_2 }^{\prime } =-\frac{g}{L}\sin \left(u_1 \right),\;\;u_2 \left(0\right)=0
% Finally, in vectorized form:
%
%
% In the script, we use odeset() to give ode45() trigger events. In this way, we can find the length of a period, by terminating the solver once u_2 (the angular velocity) falls to 0.

close all

f1 = @(u1, u2) u2; % The differential equations u' = f(t,u).
f2 = @(u1, u2, g, L) -g/L.*sin(u1);
f = @(t, u, g, L) [f1(u(1),u(2)); f2(u(1), u(2), g, L)]; % System of ode's
g = 9.81;
L = 0.1;

tspan = linspace(0, 2, 400);
phi0 = 30*pi/180; 
u0 = [phi0; 0]; % Initial conditions.
%[t U] = ode45(@(t,u)f(t,u,g,L), tspan, u0);% solve ODE for u1 and u2. This gives us angle phi and angular velocity phi'.

opts = odeset('Events',@funevent_pendulum); % Event handler. Sets ode45() to react to certain events. See funevent() at bottom of script.

phi0_vec = (10:10:170).*pi./180; % A few initial values for starting angle phi0.
T_vec = zeros(size(phi0_vec));
i = 1;
for phi0 = phi0_vec 
    [t U] = ode45(@(t,u)f(t,u,g,L), tspan, [phi0;0], opts); % Starting conditions are angle and anglular velocity.
    T_vec(i) = t(end);
    i = i + 1;

    subplot(121)
    hold on
    plot(t, U(:,1))
    xlabel('t'); ylabel('φ(t)')
    axis square
    subplot(122)
    plot(U(:,1), U(:,2))
    hold on
    xlabel('φ(t)'); ylabel('ω(t)')
    axis square
end
figure
plot(phi0_vec, T_vec)
grid on
xlabel('φ_0');
ylabel('T [s]')
title('Pendulum period T vs starting angle φ_0')

%
%% Uppgift 6: Damped pendulum
%
% The order of the diff.equation has not changed -> we still use a system of two equations.
% We simply modify the expression for {\varphi^{\prime } }^{\prime }.

close all
L = 0.1; m = 0.1; c = 0.2; g = 9.81;

f1 = @(u1, u2) u2; % The differential equations u' = f(t,u).
f2 = @(u1, u2, g, L, m) -g/L.*sin(u1) - c/m*u2;
f = @(t, u, g, L, m) [f1(u(1),u(2)); f2(u(1), u(2), g, L, m)]; % System of ode's

phi0 = deg2rad(50);
tspan = [0 1];

opts = odeset('Events',@funevent_damped_pendulum); % "Interrupt handler" for the ode solver. ie. quit when conditions met.

for phi0 = deg2rad(10:20:70)
    [t, U] = ode45(@(t,u)f(t,u,g,L,m), tspan, [phi0;0]);
    
    subplot(121)
    plot(t, U(:,1))
    hold on
    
    subplot(122)
    plot(U(:,1), U(:,2))
    hold on
    grid on
end

%
%% Randvärdesproblem/Boundary-value problems
% Here we solve 1D finite-differences boundary-value problems of the type:
%
% This is different to differential equations with only an initial value since we now have an end-bound. We solve this by forming a system of linear equations Au=b that solve u for each x.
%
% First, we divide interval [a,b] uniformly into n+1 points (n intervals). ie. we discretize the problem space.
% We approximate the 2nd derivative on the l.h.s. with a 2nd order central difference D_- D_+ u. We approximate the 1st derivative on the r.h.s with an "integer spanning" 1st order central difference D_0 u. This results in:
%
% As an example, we consider a stationary (time-invariant) heat conductor of length L along the x-axis. The temperature at each point depends on the endpoints (the boundary values).
%
% We replace the derivatives with central-difference approximations and restrict ourselves to the "inner grid" which excludes the endpoints (where the central difference is undefined).
%
% Clearly, we can now write an expression for each x_ion the inner grid.
% This results in the linear system of equations:
%
% Matrix A is sparse, in this case a tridiagonal matrix.
% We can generate it using spdiags() (sparse diagonals). Think of it as a rectangular matrix that we then shear.
% |||      \\\
% ||| =>    \\\
% |||            \\\

close all
% Problem parameters
g0 = 40;
gL = 60;
L = 1;
N = 30;
k = 2;
f = @(x) 200*exp(-(x-L/2).^2);

% Specify problem space
h = L/(N+1);
i = 0:N+1;
x = i.*h;
x_inner = x(2:end-1);

% Generate tridiagonal A coefficient matrix.
ett = ones(N,1); 
A = spdiags([-ett 2*ett -ett], [-1 0 1], N, N);
%                          ^ for row i, place the input pattern at col i-1, i, i+1
full(A); % Output full matrix (not in sparse form).

% Form b matrix, selecting only inner grid.
b = (h^2/k.*f(x_inner))';
% Insert boundary values. According to 2nd order central difference, the
% boundary values appear as elements u_0 and u_N+1 in the expressions for f_1 and f_N, respectively. We simply move them to the r.h.s. ie. matrix b. 
b(1) = b(1) + g0;
b(end) = b(end) + gL;
% Solve for u on the inner grid.
u = A\b;
% Finally complete the solution with the boundary values
u = [g0; u; gL ];
plot(x, u)
title('-ku''''=200exp(-(x-L/2)), g_0=40, g_L=60')
xlabel('x')
ylabel('Temp(x)')
grid on
axis([0 1 40 65])

%
%% 3.2 Inskjutningsmetoden/Shooting method (https://en.wikipedia.org/wiki/Shooting_method)
% In the previous boundary-value problem, we rewrote a 2nd order differential equation in terms of Euler approximations and formed a large system of equations (one for each discrete point in the problem space [excluding end-points]).
%     With this method, we reduce the boundary-value problem to an initial-value problem. We vary the initial parameters ("shooting" along different trajectories) until we find a solution that satisfies the boundary conditions. Basically, it makes systemic the process of estimating initial values, solving and adjusting until a solution is found.
%     Consider the boundary value problem from before:
% {u^{\prime } }^{\prime } \left(x\right)=f\left(x,u,u^{\prime } \right)\;\;\;a\le x\le b
% u\left(a\right)=g_a \;\;\;\;u\left(b\right)=g_b
% We can rewrite this as a system of 1st order equations:
%
% Where g() is an arbitrary function that gives 0 when the boundary conditions are satisfied.
% An equivalent initial-value problem would be:
%
% If we write the solution as a function of x and the initial value \overrightarrow{s}, \overrightarrow{u} \left(x,\overrightarrow{s} \right). That is, this is the function that will evolve as a result of the initial value.
% The solution to the boundary-value problem is then:
%
% In other words, for initial value s, the function will evolve in such a way that u\left(b,s\right)is the ending boundary, thus solving the problem.
% Numerically, we solve this by finding the root of q(s) eg. using fzero/newton's method.
%% Solve the previous heat transfer problem using the Shooting method
% A note on functions, anonymous functions and function-handles
% Anonymous function are limited to a single executable statement, beyond that we must use "real" functions.
% If "real" functions call other user defined functions, anonymous functions won't be in scope (unless passed with a function-handle), and so they must also be "real" functions.
% The '@' creates a function-handle. Anonymous functions are already defined as function-handles, but "real" functions are not, so they must be given one before they can be passed as an argument.
% "Real" functions are special in that they can be called directly from anywhere without the need for a function handle.

close all

g0 = 40;
gL = 60;
L = 1;
k = 2;

s = linspace(-10,110); % Guesses for unknown u' that will result in correct upper boundary value u(b,s)=gL. 
xspan = [0 1]; % Problem domain, NOT boundary conditions.

q = zeros(size(s));
for k = 1:length(s)
    [x, u] = ode45(@ode_system, xspan, [g0; s(k)]);
    q(k) = u(end,1) - gL; % Compare obtained boundary value with actual one.
end
plot(s,q)
hold on
title('Error in upper bound (q(s)) for various guesses for initial conditions (s) ')
xlabel('s')
ylabel('q(s)')
grid on

q_root = fzero(@(s)q_heat_transfer(s), 0) % Find initial value s that results in boundary satisfied.
plot(q_root, 0, 'ro')
figure 
[x u] = ode45(@ode_system, xspan, [g0; q_root]);
plot(x, u(:,1))
title('Heat transfer')
xlabel('x'); ylabel('u(x)')
hold on
%plot(x, u(:,2))
grid on
axis([0 1 40 65])

%

function ret = f_source(x) % functions for uppgift 8.
    ret = 200*exp(-(x-1/2).^2);
end

function q = q_heat_transfer(s)
    [x u] = ode45(@ode_system, [0 1], [40; s]);
    q = u(end,1) - 60;
end

function uprim = ode_system(x,U)
    uprim = [U(2); f_source(x)./(-2)];
end

%

function q = findq(s) % functions uppgift 9
[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

%

function [value, isterminal, direction]=funevent_pendulum(t,U)
    value = U(2); % This is the variable to react to. Angular velocity. 
    isterminal = 1; % Quit ode45() once this event occurs.
    direction = -1; % +-1: ignore unless rising/falling towards zero. 0 to allow any. We want to angular velocity to rise towards zero (it swings back in positive direction to complete the period). 
end

function [value, isterminal, direction]=funevent_damped_pendulum(t,U)
i = 4
    value = U(2); % This is the variable to react to. Angular velocity. 
    isterminal = 1; % Quit ode45() once this event occurs.
    direction = -1; % +-1: ignore unless rising/falling towards zero. 0 to allow any. We want to angular velocity to rise towards zero (it swings back in positive direction to complete the period). 
end