edolin2( f, x0, xn, y0, yn, n)

function [x, y] = edolin2( f, x0, xn, y0, yn, n)
% EDOLIN2 Rsolution d'une quation du second ordre linaire par
%         discrtisation
%         function [x, y]=edolin2( f, x0, xn, y0, yn, n)
%         l'quation est p(x)y"+q(x)y'+r(x) = q(x)
%         y(x0)=y0; y(xn)=yn;
%
%         f est le nom de la fonction qui value sur le vecteur colonne x
%         renvoie 4 colonnes formes de p(x) r(x) s(x) et q(x)
%         n est le nombre de points calculs entre x0 et xn (500 par dfaut)
if nargin <= 5
  n = 500
end;
a = zeros(n,3);
x = linspace(x0, xn, n+2);
h = (xn-x0)/(n+1);
h2 = 1./(h.*h);
hx2 = 0.5./h;
vals = feval(f, x(2:n+1));
psurh2 = vals(:,1).*h2;  %p/(h*h)
rsur2h = vals(:,2).*hx2; %r/(2h)
q      = vals(:,4);
a(:,1) = psurh2 - rsur2h;
a(:,2) = vals(:,3)- 2.0.*psurh2;
a(:,3) = psurh2 + rsur2h;
q(1) = q(1) - a(1,1)*y0;
q(n) = q(n) - a(n,3)*yn;
a(1,1)=0;
a(n,3)=0;

i = [1 2:n   ,1:n  1:n-1 3]';
j = [3 1:n-1, 1:n, 2:n   1]';
a = sparse( i, j, a(:));

if n >= 500
   y = sweep(a,q);
else
   y = a\q;
end;
y = [y0 y' yn];