% randolinj, y"+c*y^2=(6-y)x;  y(0)=1; y(4)=-1, fallet c=0.1
% bandmatrismetod, olinjärt system, utnyttja Newtons metod
  clear, clf
  xslut=4;  
  y0=1;  yslut=-1;
  xrak=[0 xslut];  yrak=[0 0];  plot(xrak,yrak); hold on
  title('y" + c*y^2 =(6-y)x,  y(0)=1, y(4)=-1 ')
  ntot=40;  
  h=xslut/ntot;  n=ntot-1; x=h*(1:n)';

% Lös först det linjära problemet med c = 0 för att skaffa
% startgissning till olinjära problemet med c > 0
  c=0
  dia=-(2-h^2*x);  
  ett=ones(n-1,1); 
  b=6*h^2*x;   
  b(1)=b(1)-y0;  b(n)=b(n)-yslut;
  y=tridia(dia,ett,ett,b);
  X=[0; x; xslut]; Y=[y0; y; yslut]; plot(X,Y,'--')
  [ymax,ix]=max(y); ymax
  text(x(ix),ymax,'c=0'), drawnow

% Nu det olinjära problemet, arbeta stegvis med ökat c eftersom 
% god startapproximation kraevs for att det ska konvergera!
  cvec=[0.025 0.05 0.1];         % experimentellt bestämda!
  for nr=1:3
    c=cvec(nr)
    relfel=1;  iter=0;
    while relfel>1e-4 & iter<15
      G=dia.*y+c*h^2*y.^2;  Gprim=dia+2*c*h^2*y;
      F=[y0; y(1:n-1)]+G+[y(2:n);yslut]-6*h^2*x;
      dy=-tridia(Gprim,ett,ett,F);       % jacobianen tridiagonal 
      relfel=norm(dy,inf)/norm(y,inf)
      y=y+dy; iter=iter+1;
    end
    if iter==15, disp('Ingen konvergens'), break, end
    Y=[y0; y; yslut];  plot(X,Y,':')
    [ymax,ix]=max(y); ymax
    text(x(ix),ymax,['c=' num2str(c)]), drawnow
  end
  plot(X,Y,'c')          % rita slutliga kurvan heldragen



