% InhomPoisson_Shooting_Example.m
%
% example code demonstrating the use of the shooting method on the Poisson
% equation
% problem: u''(x) = f, 
% a <= x <= b
% u(a) = 0, u(v) = 0
% 
% in the present test problem,
% a = 0, b = 1, f = x(x-1)
% solution: x^4/12-x^3/6+x/1
% 
% the boundary value problem is solved as an initial value problem with the
% introduction of new variables
% u1 := u
% u2 := u'
% with these, the new eigenvalue problem is
% u1'(x) = u2(x)
% u2'(x) = f(x);
% with initial condition: u1(0) = 0, u2(0) = xd;
% this can be interpreted as a nonlinear equation: the input is xd, the
% output is u1(b), and we look for a xd for which u1(b) = 0
%
% the shooting method is similar in the case of the flow stability
% equations, except in the case of flow stability, the full solution is known at one end, and alpha
% (comlex wavenumber) us tuned so that the boundary condition is satisfied
% at the other end
%%
clear all
close all
clc

%% parameters
% domain boundaries
a = 0;
b = 1;
domain = [a,b];

% options for the ODE solver: absolute and relative error tolerances
odeopt = odeset('Reltol',1e-12,'Abstol',1e-12);

% calculating source term (f) on an equidistant grid on the domain with 100 gridpoints
xvec = linspace(a,b,100).';
fvec = (xvec-1).*xvec;

%% function handles
% with these additional parameters can be passed when calling a function
%       the first one is a handle for an ordinary differential equation, which
% can be passed to an ODE solver, which is a build-in  MATLAB function for solving an Ordinary Differential Equation
%       the second one is a nonlinear equation, which is the solution of the
% Poisson equation as an initial value problem
%
% the parameters of the ODE function handle are 
%       the independent variable (x) and
%       the depedent (u(x)) variables
% here we have two additional parameters which are 'burned in' here
%    these are the grid for the source term and the source term at the gridpoints
%               for details see the function below
fhandle = @(x,u) PoissonODE(x,u,xvec,fvec);
%
%
% the parameters of the nonlinear equation function are the boundary value problem formulated as an
% ODE (fhandle), the the domain of the computation (domain)
% the only remaining input is x0
ftosolve = @(xd) Poissonshooting(fhandle,domain,xd);

% initial guess for the nonlinear problem
xd0 = 10;

% solution of the nonlinear problem - the unknown first derivative
xder = fsolve(ftosolve,xd0);

% recalculating the solution of the Poisson equation with the known proper
% initial condition
[x,sol] = ode45(fhandle,[a,b],[0;xder],odeopt);

% plotting solution
figure; plot(x,sol(:,1));
hold on;
uexact = x.^4 / 12 - x.^3 /6 + x/12;
plot(x,uexact,'-*');
xlabel('x');
ylabel('u(x)');

% plotting error
% note that the error increases linearly, as the solution error accumulates
err = abs(uexact - sol(:,1));
figure;
plot(x,err,'-*');
xlabel('x');
ylabel('abs(u_{exact}-u)');

%% functions for the nonlinear problem and the ODE

function u_b = Poissonshooting(fh,domain,xd0)
% function formulating the nonlinear problem of solving the boundary value
% problem as an initial value problem
% 
% parameters: 
% fh: function handle for the ODE solver
% domain: domain for ODE solution
% xd0: the derivative, which we are 
% 
% in order to solve with fsolve, the inputs fh and domain need to be added
% before
% this can be done by creating a function handle

% ode option object for ODE tolerances
odeopt = odeset('Reltol',1e-12,'Abstol',1e-12);

% solution of the ODE with solver ode45
% parameters: function handle, solution domain, initial condition, ode option object
% initial condition: the function is zero, the derivative is the input
[~,x_] = ode45(fh,domain,[0;xd0],odeopt);

% the output is the boundary condition at the other end: u(b) = 0
u_b = x_(end,1);
end

function dudt = PoissonODE(x,u,xvec,fvec)
% ODE function for the Poisson-equation
% u'(1) = u(2);
% u'(2) = f
% parameters:
% x: independent variable
% u: dependent variable - 2x1 column vector
% xvec: grid points for source term
% fvec: source term at grid points
% 
% in order to use with ode solver (e.g. ode45), the inputs xvec and fvec need to be added before
% this can be done by creating a function handle

dudt(2,1) = interp1(xvec,fvec,x,'pchip'); % interpolation
dudt(1,1) = u(2);
end