Changing pdepe m = 0 to m = 2 for spherical symmetry. Slight parameter adjustments.

ternary_frap.m

@@ -5,7 +5,7 @@
 
 a = -50;
 b = 0.025;
-c = 1/100000;
+c = 1/10000;
 u0 = 0.05;
 e = 0.4;
 
@@ -17,16 +17,16 @@ x = linspace(g_unique(1), g_unique(end-1), 300);
 g_inv = spacing(x, a, 2*b, e);
 g_inv = g_inv(2:end-1);
 x = [linspace(49.0, g_inv(2), 30), g_inv(3:end-2), ...
-    linspace(g_inv(end-1), 51, 30), linspace(51.1, 200, 300)];
+    linspace(g_inv(end-1), 51, 30), linspace(51.1, 300, 300)];
 %% Solve pde
 tic
 % x = [linspace(49.0, 49.8, 30), linspace(49.81, 50.2, 1000), ...
 %      linspace(50.21, 51, 30) linspace(51.01, 200, 600)];
-t = linspace(0.001, 1, 1000);
+t = linspace(0.001, 1, 100);
 fh_ic = @(x) flory_ic(x, a, u0);
 fh_bc = @(xl, ul, xr, ur, t) flory_bc(xl, ul, xr, ur, t, u0);
 fh_pde = @(x, t, u, dudx) flory_hugg_pde(x, t, u, dudx, a, b, e, c, u0);
-sol = pdepe(0, fh_pde, fh_ic, fh_bc, x, t);
+sol = pdepe(2, fh_pde, fh_ic, fh_bc, x, t);
 toc
 %% Plotting
 figure(1); hold on;
@@ -85,11 +85,11 @@ function [pl,ql,pr,qr] = flory_bc(xl, ul, xr, ur, t, u0)
 end
 
 function g0 = gamma0(x, a, b, e)
-    g0 = 10*e*(tanh((x+a)/b)+1)/2+0.001;
+    g0 = 1*e*(tanh((x+a)/b)+1)/2+0.001;
 end
 
 function sp = spacing(x, a, b, e)
-    sp = b*atanh(2/(10*e)*(x-0.001)-1)-a;
+    sp = b*atanh(2/(1*e)*(x-0.001)-1)-a;
 end
 
 function p = phi_tot(x, a, b, e, u0)