% Finite-Difference method for linear BVP.
% Solves: y''=p(t)y'+q(t)y+r(t),  y(a)=ya  y(b)=yb, t in [a,b].
% Input: p,q,r -are functions of t, [a,b] - t domain, m - number of nodes.
%        ya,yb - boundary conditions
% Output: y=(y0,...,ym) - solution.
function [t y]=fdlin(p,q,r,a,b,m,ya,yb)
h=(b-a)/m; t(1)=a; t(m+1)=b;
% Computing pk,qk,rk
for k=1:m-1
    t(k+1)=a+k*h;
    pv(k+1)=p(t(k+1));
    qv(k+1)=q(t(k+1));
    rv(k+1)=r(t(k+1));
end
% Building A and B matrices for Az=B.
A=zeros(m-1,m-1); B=zeros(m-1,1);
A(1,1)=-(4+2*h^2*qv(2)); A(1,2)=2-h*pv(2);
A(m-1,m-1)=-(4+2*h^2*qv(m)); A(m-1,m-2)=2+h*pv(m);
B(1)=2*h^2*rv(2)-(2+h*pv(2))*ya; B(m-1)=2*h^2*rv(m)-(2-h*pv(m))*yb;
for j=2:m-2
    A(j,j)=-(4+2*h^2*qv(j+1));
    A(j,j-1)=2+h*pv(j+1);
    A(j,j+1)=2-h*pv(j+1);
    B(j)=2*h^2*rv(j+1);
end
A % Display A
B % Display B
y=zeros(m+1,1);
y(1)=ya; y(m+1)=yb;
z=A\B;
for k=1:m-1;
    y(k+1)=z(k);
end
___________________________________________________________________________
% Solving y''=y+t(t-4), y(0)=y(4)=0, t in [0,4] with m=4;
clear;
p=inline('0');
q=inline('1');
r=inline('x*(x-4)');
a=0; b=4; ya=0; yb=0; m=4;
[t y]=fdlin(p,q,r,a,b,m,ya,yb)
A =

    -6     2     0
     2    -6     2
     0     2    -6


B =

    -6
    -8
    -6


t =

     0     1     2     3     4


y =

         0
    1.8571
    2.5714
    1.8571
         0