The finite element method applied to the one-dimensional example:
        ( (1+x^2) u' )' = -1 ; u(1)=u(-1)=0.
        (For simplicity we set the gravitational acceleration equal to 1.)

Begin by specifying the number of subintervals into which we divide the interval [-1,1] and the size of the subintervals.

n=20

delta=2/n

Specify the locations of the interior vertices.

Do[v[j]=-1+delta*j,{j,1,n-1}]

Define the trial functions and their derivatives.

Do[T[x_,j]=UnitStep[x-(v[j]-delta)]*UnitStep[(v[j]+delta)-x]*(1-(Abs[x-v[j]]/delta)),{j,1,n-1}]

Do[DT[x_,j]=(1/delta)*UnitStep[x-(v[j]-delta)]*UnitStep[v[j]-x]-(1/delta)*UnitStep[x-v[j]]*UnitStep[v[j]+delta-x],{j,1,n-1}]

Plot[T[x,3],{x,-1,1}]

Plot[DT[x,3],{x,-1,1}]

Define the matrix K and the vector F.

Do[F[j]=Integrate[(-1)*T[x,j],{x,-1,1}],{j,1,n-1}]

Do[Do[K[i,j]=-Integrate[DT[x,i]*DT[x,j]*(1+x^2),{x,-1,1}],{j,1,n-1}],{i,1,n-1}]

KK=Table[K[i,j],{i,n-1},{j,n-1}]

MatrixForm[KK]

FF=Table[F[j],{j,n-1}]

We now solve the matrix equation which gives us the coefficients in our approximate solution.

UU=Inverse[KK].FF

Finally, we reconstruct the approximate solution.

U[x_]=Sum[UU[[j]]*T[x,j],{j,1,n-1}];

g1=Plot[U[x],{x,-1,1}]

g2=Plot[-(1/2)*Log[(1+x^2)/2],{x,-1,1}]

Show[g1,g2]