# Solution of Cahn-Hilliard equation using Fully-Discrete scheme by Backward Euler Method

My codes:

``````mesh Th=square(100,100);
plot(Th,wait=1);
macro uex(t)(exp(-2*t)*(sin(pi*x)^2)*(sin(pi*y)^2)) //

//Fespace
fespace Vh(Th,P1);
fespace Vh2(Th,[P1,P1]);

//Parameter
real a=0.01;
real dt=0.10;
int i,k;

// Function
func real df(real u){return u^3-u; }
func real ddf(real u){return 3*u^2-1; }

Vh2 [u,w],[v,z],[oldv,oldz],[auxv,auxz],[vv,zz];
Vh dfalpha;
Vh ddfalpha;

// Problem
varf vdJ([v,z],[phi,psi])=int2d(Th)
((oldv-u)*phi+dt*(dx(oldz)*dx(phi)+dy(oldz)*dy(phi))+oldz*psi-
a*(dx(oldv)*dx(psi)+dy(oldv)*dy(psi))-dfalpha*psi);

varf vhJ([v,z],[phi,psi])=int2d(Th)
(v*phi+dt*(dx(z)*dx(phi)+dy(z)*dy(phi))+z*psi-a*(dx(v)*dx(psi)
+dy(v)*dy(psi))-ddfalpha*v*psi);

[u,w]=[(sin(pi*x)^2)*(sin(pi*y)^2),-2*a*(pi^2)*(sin(pi^2*(x^2+y^2))
-2*(pi^2)*(sin(pi*x)^2)*(sin(pi*y)^2)*(x^2+y^2))+((sin(pi*x)^2)^3)*((sin(pi*y)^2)^3)
-(sin(pi*x)^2)*(sin(pi*y)^2)];

plot(u,wait=1,cmm="Cahn-Hilliard" ,value=true);
real t;
for (int i=0;i<1./dt; i++) // Time loop
{
t=t+dt;
{
oldv[]=u[];
for (int k=0;k<10;k++)// Newton loop
{
dfalpha=df( oldv );
ddfalpha=ddf( oldv );
auxv[]=vdJ(0,Vh2);
real res= auxv[]'*auxv[]; // Residual
cout << i << " residu^2 = " << res << endl;
matrix H;
if (res< 1e-12) break;
H=vhJ(Vh2,Vh2, factorize=1,solver=LU);
vv[]=H^-1*auxv[]; // Newton
v[] -=vv[];
oldv[]=v[];
}
u[]=v[];
plot(u,wait=1,cmm="Cahn-Hilliard", value=true);
}
cout << "t" << t << "L^2-error=" << sqrt(int2d(Th)((u-uex(t))^2))<<endl;
}

**Results(or out put):**(Figures)