%{
---------------------------------------------------------------------
Solver for 1D boundary-value problem second-order ODE 
with Dirichlet boundary conditions
Method: unilizes the shooting method + the method of secants
(calls 4th-order Runge-Kutta to solve the initial value problem)
Demo version for d2y/dx2 = f(x,y,dy/dx) equation
AG: Last revision 2 April 2022
To DO later: 
1: print if no solution
2: upgrade to Neymann and mixed boundary conditions
----------------------------------------------------------------------
%}

clear  % Clears command window
clc    % Clears command history
clf    % Removes anything in the figure window before simulation
tic    % set timer for the program

f = @f1;        % function handle - pass the name of the function to shoot2
n = 51;         % number of base points
eps = 1.0e-8;  % target tolerance

x = zeros(1,n);
y = zeros(1,n);

%* print the header data
Date = today('datetime');
fprintf(' Boundary Value Problem       \n')
fprintf(' Method: Shooting + RK 4th order \n')
fprintf(' Date:   %s  \n', Date)
fprintf(' Tolerance  %4.2e \n ',eps)
fprintf('      x        y         dy \n')

% boundary values
x(1) =  0.0;
y(1) =  0.0;
x(n) =  1.0;
y(n) =  1.0;

% two assumptions for y'(1) - use dy(1) and dy(2) here only as a storage
dy(1) =  1.0;
dy(2) =  2.0;

% === call the solver
[xn,yn,dyn] = shoot2(f,x,y,dy,n,eps);
% ===

% print solutions
for i=1:n
  fprintf(' %8.4f %10.2e %10.2e \n', xn(i), yn(i), dyn(i))
end

figure (1)
plot(xn,yn,'r-o')
title('Boundary Value Problem');
xlabel('x');
ylabel('y(x)');
grid;

timeE = ceil(toc);
fprintf('\n  time elapsed = %4i sec \n',timeE)


function [f1] = f1(x,y,dy)
%{
!------------------------------------------
 the second derivative - use original ODE
 d2y/dx2 = f(x,y,dy)
------------------------------------------
%}
%f1 = -5.0*dy-4.0*y*y+exp(x);
%f1 = -4.0*dy - 6.25*y + exp(x);
f1 = 0-4.0*x*y*dy-(6.25*y*y*sin(x*y))*y+exp(x);
end


function [x,y,dy] = shoot2(f,x,y,dy,n,eps)
%{
!======================================================================
! Solution of the boundary-value second-order 1D ODE
! d2y/dx2 = f(x,y,dy/dx) with Dirichlet boundary conditions
! y(xmin) = ..., and y(xmax) = ...
! Method: unilizes the shooting method based on the method of secants
! (calls 4th-order Runge-Kutta to solve the initial value problem)
! written by:  Alex Godunov (last revision - 9 October 2009)
!----------------------------------------------------------------------
! input ...
!  f(x,y,dy)  - function d2y/dx2 (supplied by a user)
!  x(1), x(n) - boundary points
!  y(1), y(n) - boundary values (Dirichlet boundary conditions)
! dy(1),dy(2) - two guesses for y'(x(1))
!          n  - number of grid points
!        eps  - tolerance (abs value)
! output ...
!  y(i) and dy(i) solutions at points x(i) (i=1,...,n)
! note: dy corresponds to y' (the first derivative)
!======================================================================
%}

it=101;          % max number of iterations

g = zeros(1,it);
c = zeros(1,it);

%  g(it) - array of dy(it) values [iterative improvement]
%  c(it) - array of y(n) values [iterative solutions for g(it)]

% first guesses for g(it)
g(1) = dy(1);
g(2) = dy(2);

%! remember the second boundary condition 
%! since y(n) recalculated for each new g(i)
yn = y(n);

% generate base points x(i) from x(1), x(n) and n
dx = (x(n)-x(1))/(n-1);

for i=2:n
  x(i) = x(i-1)+dx;
end

% shooting iterations (for the first two - we use assumed values of dy(1))
for k=1:it
  dy(1) = g(k);
  [y,dy] = rk4_2d(f,x,y,dy,n);  % solves initial value ODE on n-base points
  c(k) = y(n);
  if abs(yn-c(k)) < eps
      break
  end
  if k >= 2
      g(k+1)=g(k)-(c(k)-yn)*(g(k)-g(k-1))/(c(k)-c(k-1));
  end
  
end

 fprintf(' Interations %3i \n',k-2)

end % function shoot2


function [y,dy] = rk4_2d(f,x,y,dy,n)
%{
================================================================
 Solution of the second-order 1D ODE (Initial-value problem)
 method:  Runge-Kutta 4th-order
 written by: Alex Godunov (last revision March 2021)
----------------------------------------------------------------
 input ...
  f(x,y,dy)- function from d2y/dx2=f(x,y,dy) (supplied by a user)
  x(i)  - array of x base point (n points)
  y(1)  - initial value for y(1)
 dy(1) - initial value for y'(1)

 output ...
  y(i)  - solutions for y  in n points
 dy(i)  - solutions for y' in n points 
=================================================
%}

for k=2:n
   h   = x(k)-x(k-1);
   k11 = h*dy(k-1);
   k12 = h*f(x(k-1),y(k-1),dy(k-1));
   k21 = h*(dy(k-1)+k12/2.0);
   k22 = h*f(x(k-1)+h/2.0, y(k-1)+k11/2.0, dy(k-1)+k12/2.0);
   k31 = h*(dy(k-1)+k22/2.0);
   k32 = h*f(x(k-1)+h/2.0, y(k-1)+k21/2.0, dy(k-1)+k22/2.0);
   k41 = h*(dy(k-1)+k32);
   k42 = h*f(x(k-1)+h,y(k-1)+k31,dy(k-1)+k32);

   y(k) = y(k-1) + (k11 + 2.0*(k21+k31) + k41)/6.0;
  dy(k) = dy(k-1)+ (k12 + 2.0*(k22+k32) + k42)/6.0; 
end
end % function rk4_2d