f1d.m

function [result,t] = fdm(a,b,h,alpha,beta,p,q,r)
    t = a+h:h:b-h;
    n = length(t);
    forward = @(t) 1 + (p(t).*(h./2));
    center  = @(t) q(t).*(h.^2) - 2;
    backward = @(t) 1 - (p(t).*(h./2));    
    right = @(t) r(t).*(h.^2);
    A = spdiags([forward(t).',center(t).',backward(t).'],[-1,0,1],n,n);
    F = right(t).';
    F(1) = F(1) - alpha;
    F(n) = F(n) - beta;
    result = A \ F; 
    result = [alpha;result;beta];
    t = [a,t,b];
end