Doesn't work also with old asymmetric meshing.

@Ternary_model/create_mesh.m

@@ -2,36 +2,32 @@ 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 = 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 = [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_temp = x_temp/2;
             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)];
+    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
     end
-% 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);
+    x = [x, x(end)+x_temp];
 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