MATLAB: Solving a system of ODE over different intervals (with different conditions on parameters)

codefunctionsindiceloopsode

Good Evening dear Matlab Community,
Now that my first part of code works, I have the ambition to make it a little complicated but I can't figure out how to implement my plan in matlab:
in the code below you can see that I'm solving a system of two differential equations. Now my challenge is that I would like to solve it with different conditions on certain parameters over the full interval l. I would like to divide the full interval l into sections i each itself divided into two zones, j.
Each section i starts with a conductive zone (j=1) of say l=3 meter for which:
Aap=lcont
Acp=Aap
hap=17
kap=0.018
And ends with a convection zone (j=2) of l=2 meter for which following is true:
Aap=2*lcont
Acp=0
hap=20
kap=0.02
After this, comes the conductive zone of the next section…
I am confused about how to define indices for the sections and zones and whether I need to include loops within the functions corresponding to my ODE in order to solve my system. I have tried (and failed) by introducing following lines within the script of my ODE functions:
for j=1 
Aap=lcont
Acp=Aap
hap=17
kap=0.018
end
for j=2
Aap=2*lcont
Acp=0
hap=20
kap=0.02
end
If I define the array J=[1 2], how can I make sure that the same indice is used to solve both ODE for each section?
Ideally I would even like to be able to build-in variation for the parameters of the conductive zone between two different sections (like for section 1 and j=1, hap=17 and for section 2 and j=1, hap=5) but I guess once I understand how to implement it for one case, I'll be able to extend it to further! 🙂
[l,y]=ode45(@myODE,[0 45],[0.48,30]);
plot(l,y(:,:))
function dy = myODE(l,y)
global u
u=y(1);
global Tp
Tp=y(2);
dTpdl = myODE1(l,u,Tp);
dudl = myODE2(l,u,Tp);
dy = [dudl;dTpdl];
end
function dudl = myODE2(l,u,Tp)
lcont=2.855;             
Psatp = 0.13*exp(18-3800/(Tp+227.03));
Phi= 1-exp(-(47*u^2+0.1*Tp*u^1.06));
Ppw= Psatp*Phi;
dudl=-(kap*Aap)*(Ppw/(Tp+273.15)-0.19);
end
function dTpdl=myODE1(l,u,Tp)
lcont=2.855;    
hcp0=0.65;                 
Psatp = 0.13*exp(18-3800/(Tp+227.03));
Phi= 1-exp(-(47*u^2+0.1*Tp*u^1.06));
Ppw= Psatp*Phi;
DHs = 0.1*(u^1.1)*((Tp+275.15)^2)*(1-Phi)/Phi;
DHvap = 1000*(2501-2.3*Tp);
DHevap = DHs+DHvap;
dTpdl=((hcp0+0.955*u)*(140-Tp)*Acp+hap*(100-Tp)*Aap+DHevap*(-kap*Aap)*(Ppw/(Tp+273)-0.19))/(1.4+4.18*u);                                                                                           
end

Best Answer

Avoid using globals to provide parameters. See <Answers: Anonymous for params>
You can integerate over different time intervals by using a loop:
tPool = [0,1,2,4,8];
y0    = 0;  % Your initial value

tAll  = [];
yAll  = [];
for k = 1:numel(tPool) - 1
    % Your parameter depending on the time step:

    switch k
        case 1
          P = 17;
        case 2
          P = 23;
        case 3
          P = 10992;
        otherwise
          P = 0;
    end
    
    % Call integrator, provide current parameter by anonymous function:

    [t,y] = ode45((t,y) @fcn(t,y,P), [tPool(k:k+1)], y0)
    
    tAll  = cat(1, tAll, t);  % Append new solution to total solution

    yAll  = cat(1, yAll, y);
    
    y0    = y(end, 1);  % Update the initial value

end
If the zones are not predefined by time steps, but by positions, use a while loop until the final time is reached and an event function set by odeset:
t0    = 0;
tEnd  = 8;
y0    = 0;  % Your initial value
tAll  = [];
yAll  = [];
while t0 < tEnd
    % Your parameter depending on the time step:
    if y(1) < 10
      P = 17;
      eventP = 20
    elseif y(1) < 20
      P = 23;
      eventP = 30;
    else   
      P = 10992;
      eventP = 100;
    end
    
    % Create new event function, which stops the integration at the
    % next step:
    opt = odeset('Events', @(t,y) @myEventsFcn(t,y,eventP));
    
    % Call integrator, provide current parameter by anonymous function:
    % tEnd is not reached, because the integration is stopped by the
    % event function except for the last step:
    [t,y] = ode45((t,y) @fcn(t,y,P), [t0, tEnd], y0, opt)
    
    tAll  = cat(1, tAll, t);  % Append new solution to total solution
    yAll  = cat(1, yAll, y);
    
    t0    = t(end);
    y0    = y(end, 1);  % Update the initial value
end
function [position,isterminal,direction] = yourEventFcn(t,y,eventP)
  position   = y(1) - eventP; % The value that we want to be zero
  isterminal = 1;   % Halt integration 
  direction  = 0;   % The zero can be approached from either direction
end
Both examples are just a demonstration and they do not run!!!
Related Question