diff --git a/ternary_frap.m b/ternary_frap.m
new file mode 100644
index 0000000000000000000000000000000000000000..3e7b8ac4e4cdbb696f93222e6be79488b45f5f80
--- /dev/null
+++ b/ternary_frap.m
@@ -0,0 +1,87 @@
+% Numerical solution of ternary FRAP model with solvent, bleached and
+% unbleached species. Model is assumed to be equilibrated
+% (bleached+unbleached=const.=pt). Then bleached species initial
+% conditions are introduced. Integration of model via pdepe.
+
+%% Solve pde
+x = linspace(0, 40, 500);
+t = linspace(0, 5000, 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);
+
+%% Plotting
+figure(1); hold on;
+for i = 1:50
+    cla; plot(x, pt(x)); plot(x, sol(i, :)); pause(0.01);
+end
+
+%% Plot and check derivatives of pt
+figure; hold on;
+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));
+
+%% 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;
+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))/...
+          (pt(x)+1);
+end
+
+function u0 = flory_ic(x)
+    if x<5
+        u0 = 0.0;
+    else
+        u0 = 0.1;
+    end
+%     u0 = pt(x)*0.5;
+end
+
+function [pl,ql,pr,qr] = flory_bc(xl,ul,xr,ur,t)
+    pl = 0;
+    ql = 1;
+    pr = 0;%ur - 0.01;
+    qr = 1;
+end
+
+function [c, f ,s] = flory_pde(x, t, u, dudx) 
+% Tentatively solve system where binary mixture is assumed, with variable
+% '1' (now called pt). Basically pt is fixed and now phi_u and
+% phi_b can swap places but nothing else.
+c = 1/pt(x);
+f = dudx;
+s = -u*lap_pt(x)/pt(x);
+end
+
+function p = pt(x)
+    p = (tanh(-(x-5)*2)+1)/2+0.1;
+end
+
+function gpt = gra_pt(x, delta)
+    gpt = (pt(x+delta)-pt(x-delta))/(2*delta);
+end
+
+function lpt = lap_pt(x, delta)
+    lpt = (gra_pt(x+delta, delta)-...
+           gra_pt(x-delta, delta))/(2*delta);
+end
+
+function glpt = gralap_pt(x, delta)
+    glpt = (lap_pt(x+delta, delta)-...
+            lap_pt(x-delta, delta))/(2*delta);
+end
+
+function llpt = laplap_pt(x, delta)
+    llpt = (gralap_pt(x+delta, delta)-...
+            gralap_pt(x-delta, delta))/(2*delta);
+end
\ No newline at end of file