How to pass solution of pdefun1 as IC-2 for pdefun2

my pde is

Given below is my effort so far...

function pdepe_philip_v3b

% 1-b Extend the above calculation to model a potential step to any chosen 
% potential assuming that the ratio of the surface concentrations of Os2+ 
% and Os3+ obeys the Nernst equation both before and after the potential
% step. Now we want to further add a potential step which will make it a
% double potential step. 

global D n F R  
n          = 1;                                  % No. of Electrons Involved
F          = 96485.3329;                         % sA/mol
R          = 8.3145;                             % kgcm^2/s^2.mol.K
D          = 5e-06;                              % cm^2/s 
Area       = 1;                                  % cm^2 
Os_bulk    = 1e-05;                              % mol/cm^3

N          = 10; 

m          = 0;                                  % Cartesian Co-ordinates

% logarithmized xmesh for semi infinite boundary condition 
x          = logspace(log10(0.000006), log10(0.06), N); x(1) = 0;


% Time Span-1
t1         = linspace(0, 10, N);                 % s      
E          = 0.8; 
sol1       = pdepe(m, @pdev3bpde, @(x)pdev3bic(x,Os_bulk), @(xl, ul, xr, ur, t)...
                   pdev3bbc(xl, ul, xr, ur, t, E, Os_bulk), x, t1);
c1         = sol1(:, :, 1); 
I_num      = ( n*F*Area) * c1(1,:) .* sqrt(D./(t1*pi));  

% I Flux For Time Span-1
figure(1); 
plot(t1, I_num, 'r--', 'LineWidth',1);
xlabel('t (s)'); ylabel('I (A/cm^2s)'); 
hold on 

% Time Span-2
t2         = linspace(10, 20, N);                % s  
c_new      = (c1(N,:))'; 
E          = -1.6; 
sol2       = pdepe(m, @pdev3bpdeNEW, @(x)pdev3bicNEW(x,c_new),@(xl, ul, xr, ur, t)...
                   pdev3bbcNEW(xl, ul, xr, ur, t, E, c_new), x, t2);

end

%% pdepe Function
%%

function [a, f, s] = pdev3bpde(x, t, u, DuDx)

global D

a           = 1; 
f           = D*DuDx;
s           = 0;

end

%% pdepe Function
%%

function [a, f, s] = pdev3bpdeNEW(x, t2, u, DuDx)

global D

a           = 1; 
f           = D*DuDx;
s           = 0;

end


%% Initial Condition-1
%%

function u0 = pdev3bic(x, Os_bulk)

u0          = Os_bulk;

end

%% Initial Condition-2
%%

function u0 = pdev3bicNEW(x, c_new)
 
u0          = c_new.*size(x);

end

%% Boundry Condition-1
%%

function [pl, ql, pr, qr] = pdev3bbc(xl, ul, xr, ur, t, E, Os_bulk)

global n F R  
 
E0          = 0.2;
alpha       = exp((E-E0).*((n*F)/(R*298.15))); 
pl          = ul - (Os_bulk./(1+alpha));
ql          = 0;
pr          = ur -  Os_bulk;
qr          = 0;     

end

%% Boundry Condition-2
%%

function [pl, ql, pr, qr] = pdev3bbcNEW(xl, ul, xr, ur, t, E, Os_bulk)

global n F R  
 
E0          = 0.2;
alpha       = exp((E-E0).*((n*F)/(R*298.15))); 
pl          = ul - (Os_bulk./(1+alpha));
ql          = 0;
pr          = ur -  Os_bulk;
qr          = 0;     

end

%%
%%

[Walter Roberson]

How does this differ from https://www.mathworks.com/matlabcentral/answers/803096-how-to-use-the-solution-of-first-pdefun-as-ic-of-the-pdefun2-while-using-pdepe?s_tid=prof_contriblnk ?

[Hashim]

It does not its just that previously i've been trying to pass the ic as a fucntion to icfun. In this attempt i've tried to make two instances of icfun one old and one "NEW". Otherwise its the same problem.

Also, I think this poster's problem is similar to mine in terms of icfun implementation but I can't seem to figure it out for the life of me.

Answer 1

[Hashim]

Just so that anybody's thinking there is a way to pass solution from pdefun1 to icfunc2... here is a solution from @Torsten on this post. I am copying the code below (it is by Torsten of course).

function main
  m = 0;
  xmesh    = linspace(0,20,101);
  tspan    = linspace(0,365*5*86400,201);
  u0       = 280.0;
  icarg    = @(x) 0.01*ones(size(x));  
  sol      = pdepe(m,@pdefun,@(x)icfun(x,icarg),@(xl,ul,xr,ur,t)bcfun(xl,ul,xr,ur,t,u0),xmesh,tspan);
  w        = sol(end,:);
  plot(xmesh,w)
  tspan2   = linspace(tspan(end),365*20*86400,201);
  u0       = 0.0;
  icarg    = @(x)interp1(xmesh,w,x);
  sol2     = pdepe(m,@pdefun,@(x)icfun(x,icarg),@(xl,ul,xr,ur,t)bcfun(xl,ul,xr,ur,t,u0),xmesh,tspan2);
  w2       = sol2(1,:);
  hold on
  plot(xmesh,w2)
end
function [c,f,s] = pdefun(xmesh,tspan,u,DuDx)
  c = 2;
  f = 1e-8*DuDx;
  s = -1e-7*u;    
end    
function u = icfun(x,icarg)
  u = icarg(x);
end
         
function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t,u0)
  pl = ul - u0;
  ql = 0;
  pr = ur;
  qr = 0;
end