%
%  lect07.m -- code for BVP example
%
%      y" + y - beta*y^3 = 0    on x = [0,L]  
%
%      y(0) = y(1) = 0
%
%  finite difference method
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  set parameters 

lambda = 0.02;  beta = 2.0;

N =  64;  L = pi*(1+lambda);  h = L/N;

%  x-axis vector & initial solution (Yj)

Xj = [h : h  : L-h]';
X  = [0 ; Xj ; L  ] ; 

Yj = sin(pi*Xj/L);

%  construct second-derivative matrix & make sparse

T  = diag(ones(N-2,1),-1) - 2*diag(ones(N-1,1),0) + diag(ones(N-2,1),1);
Ts = sparse(T);

for jj=1:1:12

%  F-functional contribution to jacobian

  Ve = (h^2)*(ones(size(Yj)) - 3*beta*(Yj.^2));

%  jacobian:  sparse solve (on/off)

  jacobn = Ts + sparse(diag(Ve,0));
%  jacobn = T  +        diag(Ve,0) ;

%  residual error:  right-side of linear solve

  residl = Ts*Yj + (h^2)*(Yj - beta*(Yj.^3));

%  solve & newton update

  flops(0);  
  deltaY = - jacobn\residl;    
  fl = flops;

  Yj = Yj + deltaY;
  disp(['  # solve ops = ' num2str(flops) ...
        '          max delta   = ' num2str(max(deltaY)) ])
end

Y = [0 ; Yj ; 0];

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
%  solution plot

clf;subplot(2,1,1)

plot(X,Y,'x')
hold on
plot(X,Y,'k')

title('\bf my first numerical nonlinear BVP')
xlabel('\bf x-axis')
ylabel('\bf y-axis')

text(1,max(Yj)*0.5,['\beta = ' num2str(beta)])
text(1,max(Yj)*0.3,['\lambda = ' num2str(lambda)])

axis([0 L 1.1*min(Y) 1.1*max(Y)])

disp(' ')
disp(['  max Yj = ' num2str(max(Yj))])
disp(['  magic constant (3/8) = ' num2str(lambda/max(Yj)^2/beta)])
disp(' ')