Analytical derivates work, parameter regime doesn't look completely off.

ternary_frap.m

@@ -4,8 +4,8 @@
 % conditions are introduced. Integration of model via pdepe.
 
 %% Solve pde
-x = linspace(0, 40, 500);
-t = linspace(0, 5000, 50);
+x = linspace(0,500, 6000);
+t = linspace(0, 500, 50);
 % opt = odeset('RelTol',1e-14, 'AbsTol', 1e-16,'MaxStep',1e-2); 
 % sol = pdepe(0,@flory_pde, @flory_ic, @flory_bc,x,t, opt);
 sol = pdepe(0,@flory_hugg_pde, @flory_ic, @flory_bc,x,t);
@@ -13,28 +13,36 @@ sol = pdepe(0,@flory_hugg_pde, @flory_ic, @flory_bc,x,t);
 %% Plotting
 figure(1); hold on;
 for i = 1:50
-    cla; plot(x, pt(x)); plot(x, sol(i, :)); pause(0.01);
+    cla; plot(x, pt(x)); plot(x, sol(i, :)); pause(0.04);
 end
 
 %% Plot and check derivatives of pt
 figure; hold on;
+x = linspace(140, 170, 500);
 plot(x, pt(x)); plot(x, gra_pt(x, 0.001)); plot(x, lap_pt(x, 0.001));
 plot(x, gralap_pt(x, 0.001)); plot(x, laplap_pt(x, 0.001));
-
+plot(x, gra_a(x, -155, 0.5));
+plot(x, lap_a(x, -155, 0.5));
+% pt(0)
+% gra_pt(0, 0.00001)
+% lap_pt(0, 0.00001)
+% lap_a(0, -155, 0.5)
+% gralap_pt(0, 0.00001)
+% laplap_pt(0, 0.00001)
 %% Function definitions for pde solver
 function [c, f ,s] = flory_hugg_pde(x, t, u, dudx)  
 % Solve with full ternary model.
-X = 2.2;
-K = 1;
+X = 3.2;
+K = 1.5;
 c = 1/(pt(x)+1);
 f = dudx;
-s = u*(lap_pt(x, 0.001)-2*X*(lap_pt(x, 0.001)-...
-        gra_pt(x, 0.001)^2-pt(x)*lap_pt(x, 0.001))+...
-        K*(-laplap_pt(x, 0.001)+...
-            gra_pt(x, 0.001)*gralap_pt(x, 0.001)+...
-            pt(x)*laplap_pt(x, 0.001)))/(pt(x)+1)+...
-    dudx*(-2*X*(gra_pt(x, 0.001)-pt(x)*gra_pt(x, 0.001))-...
-          K*gralap_pt(x, 0.001)+K*pt(x)*gralap_pt(x, 0.001))/...
+s = u*(lap_pt(x, 0.0001)-2*X*(lap_pt(x, 0.0001)-...
+        gra_pt(x, 0.0001)^2-pt(x)*lap_pt(x, 0.0001))+...
+        K*(-laplap_pt(x, 0.0001)+...
+            gra_pt(x, 0.0001)*gralap_pt(x, 0.0001)+...
+            pt(x)*laplap_pt(x, 0.0001)))/(pt(x)+1)+...
+    dudx*(-2*X*(gra_pt(x, 0.0001)-pt(x)*gra_pt(x, 0.0001))-...
+          K*gralap_pt(x, 0.0001)+K*pt(x)*gralap_pt(x, 0.0001))/...
           (pt(x)+1);
 end
 
@@ -42,7 +50,7 @@ function u0 = flory_ic(x)
     if x<5
         u0 = 0.0;
     else
-        u0 = 0.1;
+        u0 = 0.5;
     end
 %     u0 = pt(x)*0.5;
 end
@@ -64,7 +72,7 @@ s = -u*lap_pt(x)/pt(x);
 end
 
 function p = pt(x)
-    p = (tanh(-(x-5)*2)+1)/2+0.1;
+    p = (tanh(-(x-155)*2)+1)/2+0.3;
 end
 
 function gpt = gra_pt(x, delta)
@@ -84,4 +92,12 @@ end
 function llpt = laplap_pt(x, delta)
     llpt = (gralap_pt(x+delta, delta)-...
             gralap_pt(x-delta, delta))/(2*delta);
+end
+
+function grad = gra_a(x, a, b)
+    grad = -(1-tanh(-(x+a)/b).^2)*1/b*0.5;
+end
+
+function lap = lap_a(x, a, b)
+    lap = -2*tanh(-(x+a)/b).*(1-tanh(-(x+a)/b).^2)*1/b^2*0.5;
 end
\ No newline at end of file