
Clear[phi];
phi[0,x_,{t0_,t1_}]:=(x-t1)/(t0-t1);
phi[1,x_,{t0_,t1_}]:=(x-t0)/(t1-t0);

Clear[basisVec];
basisVec::usage = "basisVec[grid,x] gibt die low-order Basisfunktionen auf grid zurueck; wird verwendet um die Loesung plotten zu koennen"
basisVec[grid_,x_]:=Join[{Piecewise[{{phi[0,x,{grid[[1]],grid[[2]]}],grid[[1]]<=x<grid[[2]]}}]}, Table[Piecewise[{{phi[1,x,{grid[[k-1]],grid[[k]]}],grid[[k-1]]<=x<grid[[k]]},{phi[0,x,{grid[[k]],grid[[k+1]]}],grid[[k]]<=x<grid[[k+1]]}}],{k,2,Length[grid]-1}], {Piecewise[{{phi[1,x,{grid[[Length[grid]-1]],Last[grid]}],grid[[Length[grid]-1]]<=x<=Last[grid]}}]}];

Clear[MassLocal];
MassLocal[{t0_,t1_}]:=(t1-t0)/6 {{2,1},{1,2}};

Clear[StiffnessLocal];
StiffnessLocal[{t0_,t1_}]:= 1/(t1-t0) {{1,-1},{-1,1}};

Clear[MassGlobal];
MassGlobal[grid_]:=Module[{dim = Length[grid], mass, localmass, i, j},

   mass = Table[0,{dim},{dim}];

   For[i=1,i<=dim-1,i++,
         localmass = MassLocal[{grid[[i]],grid[[i+1]]}];
         mass[[i,i]] += localmass[[1,1]];
         mass[[i+1,i]] += localmass[[2,1]];
         mass[[i,i+1]] += localmass[[1,2]];
         mass[[i+1,i+1]] += localmass[[2,2]];
   ];
   Return[mass];
];

Clear[StiffnessGlobal];
StiffnessGlobal[grid_]:=Module[{dim = Length[grid], stiff, localstiff, i},

   stiff = Table[0,{dim},{dim}];

   For[i=1,i<=dim-1,i++,
         localstiff = StiffnessLocal[{grid[[i]],grid[[i+1]]}];
         stiff[[i,i]] += localstiff[[1,1]];
         stiff[[i+1,i]] += localstiff[[2,1]];
         stiff[[i,i+1]] += localstiff[[1,2]];
         stiff[[i+1,i+1]] += localstiff[[2,2]];
   ];
   Return[stiff];
];

Clear[LocalVec];
LocalVec[f_,{t0_,t1_}]:=Table[NIntegrate[phi[i1,x,{t0,t1}] f[x],{x,t0,t1}],{i1,0,1}];

Clear[GlobalVec];
GlobalVec[f_,grid_]:=Module[{dim = Length[grid], rhs, localrhs, i},

   rhs = Table[0,{dim}];

   For[i=1,i<=dim-1,i++,
      localrhs = LocalVec[f,{grid[[i]],grid[[i+1]]}];
      rhs[[i]] += localrhs[[1]];
      rhs[[i+1]] += localrhs[[2]];
   ];

   Return[rhs];
];

Clear[FEMSolve];
FEMSolve::usage="FEMSolve[{f,mu,lam},grid,x,{u0,u1}] berechnet FEM Loesung der ODE -mu u''[x] + lam u[x] == f[x], u[0]==u0, u[1]==u1"
FEMSolve[{f_,mu_,lam_},grid_,x_,{u0_,u1_}]:=Module[{uvec,Amat,rhs,dim = Length[grid],sol,g},

   Amat = mu StiffnessGlobal[grid] + lam MassGlobal[grid];

   If [u0!=0 || u1!=0,
      g = (f[#]-lam u0 phi[0,#,{First[grid],Last[grid]}]-lam u1 phi[1,#,{First[grid],Last[grid]}])&,
      g = f
      ];

   rhs = GlobalVec[g,grid];


   Amat = Take[Amat,{2,dim-1},{2,dim-1}];

   rhs = Take[rhs,{2,dim-1}];
   
   uvec = LinearSolve[Amat,rhs];

   sol = uvec.Take[basisVec[grid,x],{2,dim-1}];

   If [u0!=0 || u1!=0,
      sol += u0 phi[0,x,{First[grid],Last[grid]}] + u1 phi[1,x,{First[grid],Last[grid]}];
      ];

   Return[sol];

];