{   prog7a.pas -- Solve the two-point boundary value problem

        y''+P(x)y'+Q(x)y=R(x), y(x0)=y0, y(xn)=yn

    by the central difference method.
    Written May 1, 2001 by Eric Olson for Mathematics 483

    For details see Section 5.2 in Numerical Analysis for Applied
    Mathematics, Science and Engineering by Greenspan and Casulli. }

Program main;

const Nmax=1024;
type vector=array[1..Nmax] of real;
    tridiag=array[1..3] of vector;

function P(x:real):real;
begin
    P:=cos(x);
end;

function Q(x:real):real;
begin
    Q:=-(x*x+1);
end;

function R(x:real):real;
begin
    R:=2;
end;

procedure trisolve(var a:tridiag;var b:vector;n:integer;var y:vector);
    var i:integer;
        p:real;
begin
    for i:=1 to n-2 do begin
        p:=a[1][i+1]/a[2][i];
        a[2][i+1]:=a[2][i+1]-p*a[3][i];
        b[i+1]:=b[i+1]-p*b[i];
    end;
    for i:=n-1 downto 2 do begin
        p:=a[3][i-1]/a[2][i];
        b[i-1]:=b[i-1]-p*b[i];
    end;
    for i:=1 to n-1 do begin
        y[i]:=b[i]/a[2][i];
    end;
end;

procedure CD(x0,y0,xn,yn:real;n:integer;var y:vector);
    var a:tridiag;
        b:vector;
        x,h,hh2:real;
        i:integer;
begin
    h:=(xn-x0)/n;
    hh2:=2.0*h*h;
    for i:=1 to n-1 do begin
        x:=x0+i*h;
        b[i]:=hh2*R(x);
        a[1][i]:=2-h*P(x);
        a[2][i]:=-4.0+hh2*Q(x);
        a[3][i]:=2+h*P(x);
    end;
    b[1]:=b[1]-a[1][1]*y0;
    b[n-1]:=b[n-1]-a[3][n-1]*yn;

    trisolve(a,b,n,y);

    writeln('#xi yi');
    writeln(x0,' ',y0);
    for i:=1 to n-1 do begin
        x:=x0+i*h;
        writeln(x,' ',y[i]);
    end;
    writeln(xn,' ',yn);
    writeln;
end;

var i,n:integer;
    y:vector;
    x0,y0,xn,yn:real;
begin
    x0:=-2.0; y0:=-1.0;
    xn:=2.0; yn:=-1.0;

    n:=4;
    for i:=2 to 8 do begin
        CD(x0,y0,xn,yn,n,y);
        n:=n*2;
    end;
end.