Consideramos o problema
|  |  |  |  | (1.163) | 
|  |  |  |  | (1.164) | 
Aqui, computamos aproximações de elementos finitos no espaço das funções lineares por partes  com sucessivos refinamentos de malha. Utilizamos uma malha inicial uniforme com 10 células e fazemos, então, 5 refinamentos sucessivos utilizando como critério de refinamento a estratégia (1.162) com . A Figura 1.7 apresenta o esboço do gráfico da solução de elementos finitos na malha mais refinada. Além disso, na Tabela 1.1 temos os o número de células e o  máximo respectivo.
 
Com o FEniCS, a computação do problema de elementos finitos pode ser feita com o seguinte código:
from __future__ import print_function, division
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
# malha
mesh = IntervalMesh(10,0,1)
# espaco
V = FunctionSpace(mesh, ’P’, 1)
# fonte
f = Expression(’exp(-100*pow(fabs(x[0]-0.5),2))’,degree=1)
# condicoes de contorno
def boundary(x,on_boundary):
    return on_boundary
#iteracoes
for iter in np.arange(6):
    #problema
    bc = DirichletBC(V,Constant(0.0),boundary)
    u = TrialFunction(V)
    v = TestFunction(V)
    a = u.dx(0)*v.dx(0)*dx
    L = f*v*dx
    #resolve
    u = Function(V)
    solve(a == L, u, bc)
    #grafico
    plt.close(’all’)
    xx = mesh.coordinates()[:,0]
    sorted_indices =  np.argsort(xx)
    yy = u.compute_vertex_values()
    plt.plot(xx[sorted_indices],yy[sorted_indices],
                 marker="o",label=r"$u_h$")
    plt.legend(numpoints=1)
    plt.grid(’on’)
    plt.show()
    DG = FunctionSpace(mesh, "DG", 0)
    v = TestFunction(DG)
    a = CellVolume(mesh)
    eta = assemble(f**2*v*a*dx)
    # refinamento da malha
    cell_markers = MeshFunction("bool", mesh, mesh.topology().dim(), False)
    eta_max = np.amax(eta[:])
    print(eta_max)
    print("%d %d %1.1E\n" % (iter,mesh.num_cells(),eta_max))
    alpha = 0.5
    for i,cell in enumerate(cells(mesh)):
        if (eta[i] > alpha*eta_max):
            cell_markers[cell] = True
    mesh = refine(mesh, cell_markers)
    V = FunctionSpace(mesh, ’P’, 1)