lambda=D*tau0/(0.6267*W0*W0);
[phi,T] = nucleus(Nx,Ny,seed);
txx(i,j) = (T(ip,j) -2*T(i,j)+ T(im,j))./dx^2;
tyy(i,j) = (T(i,jp) -2*T(i,j)+ T(i,jm))./dy^2;
lap_T(i,j)= txx(i,j)+tyy(i,j);
nx(i,j) = (phi(ip,j) - phi(im,j))/(2*dx);
ny(i,j) = (phi(i,jp) - phi(i,jm))/(2*dy);
nxy(i,j) = (phi(ip,jp) - phi(ip,jm)- phi(im,jp)+ phi(im,jm))/(4*dy*dx);
nxx(i,j) = (phi(ip,j) -2*phi(i,j)+ phi(im,j))/(dx*dx);
nyy(i,j) = (phi(i,jp) -2*phi(i,j)+ phi(i,jm))/(dy*dy);
theta(i,j)=atan2(ny(i,j),nx(i,j));
a(i,j) = 1.0 + epsilon_m*cos(m*(theta(i,j)-theta0));
tau(i,j)= tau0*a(i,j)*a(i,j);
W_dthehta(i,j) = -W0*m*epsilon_m*sin(m*(theta(i,j)-theta0));
term0(i,j)= dt./tau(i,j);
term1(i,j)= (phi(i,j)-lambda.*T(i,j)*(1-(phi(i,j)*phi(i,j))))*(1-(phi(i,j)*phi(i,j)));
term2(i,j)= W_dthehta(i,j)*W_dthehta(i,j)+W(i,j)*m*m*(W0-W(i,j));
term3(i,j)= (2*nx(i,j)*ny(i,j)*nxy(i,j)-(nx(i,j)*nx(i,j)*nyy(i,j)+ny(i,j)*ny(i,j)*nxx(i,j)))/(nx(i,j)*nx(i,j)+ny(i,j)*ny(i,j));
term4(i,j)= 2*W(i,j)*W_dthehta(i,j);
term5(i,j)=(nxy(i,j)*(nx(i,j)*nx(i,j)-ny(i,j)*ny(i,j))-nx(i,j)*ny(i,j)*(nxx(i,j)-nyy(i,j)))/(nx(i,j)*nx(i,j)+ny(i,j)*ny(i,j));
term6(i,j)=W(i,j)*W(i,j);
term7(i,j)=nxx(i,j)+nyy(i,j);
phi(i,j)=phi(i,j)+term0(i,j).*(term1(i,j)+term2(i,j).*term3(i,j)+term4(i,j).*term5(i,j)+term6(i,j).*term7(i,j));
T(i,j)=T(i,j)+dt.*(D.*lap_T(i,j)+0.5.*(phi(i,j)-phiold(i,j)));
if(mod(istep,nprint)== 0)
fprintf('donestep:%5d\n',istep);
write_vtk_grid_values(Nx,Ny,dx,dy,istep,phi,T);
compute_time=etime(clock(),time0);
fprintf('Compute Time: %10d\n',compute_time);