%%
% Sauer Computer Problem 7.2.3 (a)
%
% Lös randvärdesproblemet y'' = 18*y^2,
%                         y(1) = 1/3, y(2) = 1/12
% med finita differensmetoden.
%
% Joar Bagge 2016-11-29

function Sauer_7_2_3_a
  % Randvillkor
  xL = 1; yL = 1/3;  % vänster rand
  xR = 2; yR = 1/12; % höger rand

  % Dela upp intervallet i n delintervall
  figure(1), clf
  figure(2), clf
  for n=[10, 20, 40];
    x = linspace(xL, xR, n+1);
    h = (xR-xL)/n; % delintervallens längd

    % Approximera derivatorna med centraldifferenser...
    % detta leder till ett ickelinjärt system f(y)=0.

    % Lös med Newtons metod
    y = zeros(n+1, 1); % startgissning
    dy = Inf;
    while norm(dy) >= 1e-12
      dy = -J(y, h)\f(y, yL, yR, h);
      y = y + dy;
    end

    % Rita lösningen
    figure(1)
    plot(x, y, '*-')
    hold on
    title('Lösning')

    % Rita felet som funktion av x
    yexact = 1./(3*x.^2)';
    abserr = abs(y - yexact);
    figure(2)
    semilogy(x, abserr, '*-')
    hold on
    title('Fel')
    MaximaltFel = max(abserr)

    pause
  end
end

function fy = f(y, yL, yR, h)
% Sätt upp vektorn f som funktion av y.
  n = numel(y) - 1;
  fy = zeros(n+1, 1);
  % Första och sista elementet:
  fy(1) = y(1) - yL;
  fy(end) = y(end) - yR;
  % De andra elementen:
  for i=2:n
    fy(i) = y(i+1) - 2*y(i) + y(i-1) - 18*h^2*y(i)^2;
  end
end

function Jy = J(y, h)
% Sätt upp Jacobimatrisen J som funktion av y.
  n = numel(y) - 1;
  Jy = zeros(n+1, n+1);
  % Första och sista raden:
  Jy(1,1) = 1; Jy(end,end) = 1;
  % De andra raderna:
  for i=2:n
    Jy(i,i-1) = 1;
    Jy(i,i) = -2 - 18*h^2 * 2*y(i);
    Jy(i,i+1) = 1;
  end
end
