%% Useful built-in Matlab methods
%
% Asgeir Birkisson, October 2015
%
% Matlab offers a plenty of useful methods for solving various problems in
% scientific computing -- quadrature (integration), root-finding, IVPs, BVPs.
% Here are examples of some of them

%% Quadrature
%
% The term for numerical intergration is quadrature. Matlab offers a
% selection of methods for quadrature, e.g. quad, quadgk, and integral. We pass
% the method an anonymous function we wish to integrate, and an interval
% (tolerance optional):
f = @(x) x.^2;
intval = integral(f, 0, 2)
%%
% Compare with exact answers
intval - (1/3)*2^3

%% Root-finding
%
% We use the method fzero to find a root of a scalar-valued function f:
f = @(x) cos(exp(x)) + exp(cos(x))
fplot(f)
%%
r = fzero(f,4)
f(r)
hold on, plot(r,f(r),'r.','MarkerSize',10)
%% Initial value problems
%
% We can use the methods ODE45 or ODE115s (for stiff problems) to solve IVPs.
% Need to rewrite as a first order system, so
%    y'' + (y^2-1)y' + y = 0
% becomes
dydt = @(t,y) [y(2); (1-y(1)^2)*y(2)-y(1)];
%%
% Need to pass the interval we want to obtain the solution on, and a vector of
% initial conditions for y and y':
y0 = 2;
yp0 = 0;
timeInterval = [0 20];
[t,y]=ode45(dydt,timeInterval,[y0 yp0]);
%%
plot(t,y,'linewidth',2);shg
title('Solution of an initial value problem')
%%
sol = ode113(dydt,timeInterval,[y0 yp0])
deval(sol, 4.4:0.01:4.5)
%% Boundary value problems
%
% Can use methods BVP4C and BVP5C for solving boundary value problems. Again,
% need to rewrite as a first order system, so
%   y'' + sin(y) = 0
% becomes
dydt = @(t,y) [y(2); -sin(y(1))];
%%
% For boundary conditions, also need to in separate components
bcs = @(ya,yb) [ ya(1)-2; yb(2)+1];
%%
% Need to pass initial guesses through the bvpinit method
solinit = bvpinit(linspace(0,4,50), [1 0])
%%
% The solution will be a Matlab struct
sol = bvp5c(dydt, bcs, solinit)
%%
% We can for example look at what gridpoints the solution has been computed at:
sol.x
%%
% Plot the solution
plot(sol.x, sol.y,'linewidth',2)
title('Solution of a boundary-value problem'), shg
%%
% See also deval
xx = linspace(0,4,200);
yy = deval(sol, xx);
plot(xx, yy(1,:),'linewidth',2)
title('Solution of a boundary-value problem'), shg