Reverting create_mesh to symmetric.

@Ternary_model/create_mesh.m

@@ -2,32 +2,36 @@ function create_mesh(T)
 %CREATE_MESH creates mesh for solving Ternary_model.
 %   Mesh density larger for steeper gradients.
 
+if strcmp(T.mode, 'Constituent')
 % Start with left side of mesh, at r=0/x=0, with a given step size
-x = [0, 0.01/T.precision];%/10
-x = [0, 0.01/T.precision/10];%/10
-
-while x<-2*T.a
-    x_temp = x(end)-x(end-1);
-    phi_nm1 = T.phi_tot(x(end), T.a, T.b, T.e, T.u0, 0);
-    phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0);
-    if x(end) <= -T.a
-        while phi_nm1-phi_n > 0.0002/T.precision
+    x = linspace(0, 0.75, 40);
+    x = [x, x(end)+0.001/T.precision];
+    while x(end)<-T.a
+        x_temp = x(end)-x(end-1);
+        phi_nm1 = T.phi_tot(x(end), T.a, T.b, T.e, T.u0, 0);
+        phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0);
+        while phi_nm1-phi_n > 0.002/T.precision
             x_temp = x_temp/2;
             phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0);
         end
-    elseif x(end) > -T.a
-        while (phi_nm1-phi_n < 0.0001/T.precision) && ...
-                (x_temp < 0.01/T.precision)
-            x_temp = 2*x_temp;
-            phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0);
-        end
+        x = [x, x(end)+x_temp];
+    end
+    x = [x, x(end)+cumsum(fliplr(diff(x)))];
+    x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 100*T.precision);
+    x_right = linspace(-10*T.a+0.1, T.system_size, 100*T.precision);
+    T.x = [x, x_middle, x_right];
+    T.x = T.x(1:find(T.x-T.system_size>=0, 1)-1);
+
+    % Fine mesh close to x0 if ic = 'Gauss'
+    [~, ind] = min(abs(T.x-T.x0));
+    if (T.x(ind+20)-T.x(ind) > 0.1) && strcmp(T.ic, 'Gauss')
+        T.x = [T.x(1:ind-21), linspace(T.x(ind-20), T.x(ind+20),...
+                                       ceil((T.x(ind+20)-T.x(ind-20))*200)),...
+               T.x(ind+21:end)];
     end
-    x = [x, x(end)+x_temp];
+% T.x = [linspace(0, 0.99, 100), linspace(0.9901, 1.01, 2000), linspace(1.0161, 2, 100),...
+%        linspace(2.01, 10, 100), linspace(10.01, 300, 200)];
+elseif strcmp(T.mode, 'Client')
+    T.x = linspace(0, 20, 24000);
 end
-x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 200*T.precision);
-x_right = linspace(-10*T.a+0.1, T.system_size, 200*T.precision);
-x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 600*T.precision);
-x_right = linspace(-10*T.a+0.1, T.system_size, 600*T.precision);
-T.x = [x, x_middle, x_right];
-T.x = T.x(1:find(T.x-T.system_size>=0, 1)-1);
 end