BeginPackage["BVP`"];
BVP::usage="BVP[a,b,f,xvalues,BCType,BCValues] 
solves a second order boundary value problem";

Begin["`Private`"];
M0={{0.6872983346207417, 0.4,  -0.08729833462074164},  {0, 1, 0}, 
  {-0.08729833462074164,  0.4, 0.6872983346207417}};
M1={{-1.274596669241483, 1.549193338482966, -0.2745966692414834}, {-(1/2), 0, 1/2},
   {0.2745966692414834, -1.549193338482966, 1.274596669241483}};

ElementStiffnessMatrix[a_,b_,{xl_,xr_}] :=
  Module[{pl,pc,pr,h,avalues,bvalues},
     pc = (xl+xr)/2;
     h=(xr-xl)/2;     
     {pl,pr}={pc-Sqrt[3/5]*h,pc+Sqrt[3/5]*h};
     avalues=Map[a,{pl,pc,pr}];
     bvalues=Map[b,{pl,pc,pr}];
     Return[Transpose[M1].DiagonalMatrix[avalues*{5,8,5}/9].M1/h+
            Transpose[M0].DiagonalMatrix[bvalues*{5,8,5}/9*h].M0];
]
ElementRHS[f_,{xl_,xr_}] :=
   Module[{pl,pc,pr,w,h,fvalues},
     pc = (xl+xr)/2;
     h=(xr-xl)/2;     
     {pl,pr}={pc-Sqrt[3/5]*h,pc+Sqrt[3/5]*h};
     fvalues=Map[f,{pl,pc,pr}];
     Return[-Transpose[M0].(fvalues*{5,8,5})/9*h]
]
BVP[a_,b_,f_,xvalues_,BCType_,BCValues_]:=
    Module[{n,x,u,M,mElement,RHS,mRHS},
      n=Length[xvalues];
      x=Table[0,{2*n-1}];
      For[e=1,e<=n,e++,x[[2*e-1]]=xvalues[[e]]];
         For[e=1,e<n,e++, x[[2*e]]=(xvalues[[e]]+xvalues[[e+1]])/2];
      M=Table[0,{2*n-1},{2*n-1}];       
      RHS=Table[0,{2*n-1}];
      For[e=1,e<n,e++,
        mElement=N[ElementStiffnessMatrix[a,b,{xvalues[[e]],xvalues[[e+1]]}]];
        For[i=1,i<=3,i++,
          For[j=1,j<=3,j++, M[[(e-1)*2+i,(e-1)*2+j]]+=mElement[[i,j]]]];
        mRHS=N[ElementRHS[f,{xvalues[[e]],xvalues[[e+1]]}]];
        For[i=1,i<=3,i++, RHS[[(e-1)*2+i]]+=mRHS[[i]]];
      ];
      Which[BCType[[1]]=="D"||BCType[[1]]=="d",
        M[[1]]=Table[0,{2*n-1}];M[[1,1]]=1;RHS[[1]]=BCValues[[1]],
        BCType[[1]]=="N"||BCType[[1]]=="n", RHS[[1]]-=BCValues[[1]]];
      Which[BCType[[2]]=="D"||BCType[[2]]=="d",
        M[[2*n-1]]=Table[0,{2*n-1}];M[[2*n-1,2*n-1]]=1;RHS[[2*n-1]]=BCValues[[2]],
        BCType[[2]]=="N"||BCType[[2]]=="n", RHS[[2*n-1]]-=BCValues[[2]]];
      u=LinearSolve[M,RHS];
      Return[Transpose[{x,u}]];
];
End[];
EndPackage[];
