function [x u]=BVPsolver(N,F,G,a,b,l)
    //This Scilab function requires inputs [p(x), q(x), r(x), s(x)]=F(x); 
    //alpha0=G(1); alpha1=G(2); alpha2=G(3);
    //beta0=G(4); beta1=G(5); beta4=G(6);
    //gamma0=G(7); gamma1=G(8); gamma2=G(9);
    //delta0=G(10); delta1=G(11); delta2=G(12);
    //
    //Where the equation that's being solved is:
    //
    //p(x) d^2 u/dx^2+q(x) du/dx+r(x)u=s(x)
    //
    //cond1: alpha0*u(a)+alpha1*u'(a)=alpha2; 
    //
    //cond2: beta0*u(a)+beta1*u'(a)=beta2;
    //
    //cond3: gamma0*u(b)+gamma1*u'(b)=gamma2;
    //
    //cond4: delta0*u(b)+delta1*u'(b)=delta2;
    //
    //Either cond1 and cond 4 are imposed or cond 1 and 2 are imposed or cond 3 and 4 are imposed
    //
    //N+1 is the amount of Chebyshev Extrema Grid points used in the approximation
    //
    //l is a parameter such that if it equals any value but 1 will cause x to be equally spaced and u be u(x) on those equally spaced points
    //if l equals 1 then x is the chebyshev extrema grid and u is u(x) at the chebyshev extrema grid
    //For this function to work one must previously have downloaded fct2, ifct and Pseudo files from the file exchange
    
    
    [p, q, r, s]=F(x);
    alpha0=G(1); alpha1=G(2); alpha2=G(3);
    beta0=G(4); beta1=G(5); beta4=G(6);
    [ddT, dT, T, x]=Pseudo(N,a,b);
    K=diag(p)*ddT+diag(q)*dT+diag(r)*T; L=s;
    
    H=fct(H); G=fct(L);
    H(N,:)=alpha1*dT(1,:)+gamma1*dT(1,:); H(N,1)=H(N,1)+alpha0+gamma0;
    
    H(N+1,:)=beta1*dT(1,:)+delta1*dT(N+1,:); H(N+1,N+1)=H(N+1,N+1)+beta0+delta0;
    
    G(N)=alpha2+gamma2; G(N+1)=beta2+delta2;
    aa=H\G;
    if l=1 then
        u=ifct(aa);
    else
        xx=linspace(-1,1,N+1)';
        x=(b-a)/2*xx+(a+b)/2;
        TT=cos(acos(xx)*(0:N));
        u=TT*aa;
    end
endfunction