Sample Navier-Stokes Solver

Dear Friends,
I want to write a sample Navier Stokes solver using Rubin-Graves type linearization technique as u^(n+1)=u^n+u. The following is my code. I have some difficulties in the starting boundary values and also in the linearized non-linear terms. Would you please share your ideas and solutions.
Best regards,
Selcuk Han AYDIN, PhD.
===================== The Code ==================
// Parameters
int nn = 4;
real ut=1e-30; //If I set ut=1 results becomoes 1+e30
// Mesh
mesh Th = square(nn, nn);
// Fespace
// fespace Uh(Th, P1b);
fespace Uh(Th, P2);
Uh ux, uy, vx, vy, uxn, uyn;
fespace Ph(Th, P1);
Ph p, q, pn;
uxn=0; uyn=0; pn=0;

varf prob0(uxn,vn) = on(3,uxn=ut);
real[int] res=prob0(0,Uh);
cout << uxn<<“:”<<uyn<< endl;

macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [ [ux,uy]'[dx(vx),dy(vx)] , [ux,uy]'[dx(vy),dy(vy)] ]//

problem LinNS([ux,uy,p],[vx,vy,q]) =
+ UgradV(uxn,uyn, ux, uy)'
[vx,vy] //Error in these terms
+ UgradV(ux,uy,uxn,uyn)’
+ UgradV(uxn,uyn,uxn,uyn)’
- (p+pn) * Divergence(vx, vy)
+ Divergence(ux,uy)*q + Divergence(uxn,uyn)q
(p+pn)*q )
+on(1,2,3,4, ux=0., uy=0.);

What are the reasons of the errors in the following code.

  1. Why there are -nan in the first case
  2. Why there are e+29 in the second case

Best regards;

mesh Th = square(4, 4);
fespace Vh(Th, P1);
Vh u, v;
func ue = (x^2+y^2-1)/4; //ue: exact solution

solve p0(u,v) = on(1,2,3,4,u=ue)+on(0,u=0);

cout << “Solution u=” << u << endl;

varf prob0(u,v) = on(1,2,3,4,u=ue);
real[int] res=prob0(0,Vh);
cout << “u=”<<u<< endl;