1 real c = 6.25;
 2 int N  = 80;
  
 4 mesh M;
 5 {
 6     border a(t=0,1) { x=t;   y=0;   label=1; };
 7     border b(t=0,1) { x=1;   y=t;   label=1; };
 8     border c(t=0,1) { x=1-t; y=1;   label=1; };
 9     border d(t=0,1) { x=0;   y=1-t; label=1; };
10     M = buildmesh(a(N)+b(N)+c(N)+d(N));
11 }
 
13 load "Element_P3"
14 func Pk = P3;
15 fespace Vh(M, Pk);
16 Vh u;
17 func BC = cos(pi*x)*cos(pi*y);
18 varf vInit(w, v) = on(1, w = BC);
19 varf vJ(w, v) = int2d(M)(dx(w)*dx(v) + dy(w)*dy(v) - c*exp(u)*w*v) + on(1, w = 0);
20 varf vRes(w, v) = int2d(M)(dx(u)*dx(v) + dy(u)*dy(v) - c*exp(u)*v) + on(1, w = u);
21 int it = 20;
22 real tol = 1.0e-6;
23 u[] = vInit(0, Vh, tgv = -1);
24 real[int] b = u[];
25 for(int i = 1; i < it; ++i) {
26     real[int] v = vRes(0, Vh, tgv = -1);
27     v -= b;
28     real res = sqrt(v' * v);
29     cout << "Norm of the residual = " << res << endl;
30     if(res < tol)
31         break;
32     matrix J = vJ(Vh, Vh, tgv = -1);
33     real[int] inc = J^-1 * v;
34     u[] -= inc;
35     plot(u, cmm = "Solution at iteration #" + i);
36 }
37 plot(u, cmm = "Final solution");