Code for boundary value problem

I'm trying to answer this question:

This is the code I have produced:

    function [x,y] = bvp(coeff, bcs, xspan, h)
    N = (xspan(2)-xspan(1))/(h); 
    x = xspan(1):h:xspan(2);

    main = [ones(N-1,1).*(coeff(3)*h^2-2*coeff(1))];
    sub = ones(N-2,1).*(coeff(1)-coeff(2)*h/2);
    super = ones(N-2,1).*(coeff(1)+coeff(2)*h/2); 

    A = (diag(main) + diag(sub,-1) + diag(super,1));
    B = [(-(coeff(1)-coeff(2)*h/2)*bcs(1));zeros(N-3,1);-(coeff(1)+coeff(2)*h/2)*bcs(2)];
    y = A\B;
    y = [bcs(1); y; bcs(2)];
end

I'm not sure why this is not working. Any help would be appreciated.