Solving a system of ODEs whose coefficients are piecewise functions

I try to plot the solution of a system of ODE, on [-10,10], for the initial data [0.001 0.001], using the function:

function dwdt=systode(t,w)
if 0< t<1 
    f = t*(3-2*t);
    if -1<t< 0
        f=t*(3+2*t);
    else
        f = 1/t;
    end;
    if 0< t <1 
        h=4*t^4-12*t^3+9*t^2-4*t+3;
        if -1< t < 0
            h=4*t^4+12*t^3+9*t^2+4*t+3;
        else 
            h=0;
        end;  
        beta=0.5+exp(-abs(t));
        dwdt=zeros(2,1);
        dwdt(1)=-f*w(1)+w(2);
        dwdt(2)=-beta*w(1)-f*w(2)+h*w(1)-f*w(1)^2;
    end 

The coefficients f(t) and g(t) are piecewise functions as follows.

With the commands

    tspan = [-10 10];
    z0=[0.001 0.001];
    [t,z] = ode45(@(t,z) systode(t,z), tspan, z0);
    figure
    plot(t,z(:,1),'r');

I get the message

tspan = [-10 10];
 ↑
Error: Invalid use of operator. 

Where could be the mistake? I am also not sure that I defined correctly the functions f, h, β.

Answer 1

[Torsten]

function main
  T = [];
  Z = [];
  z0=[0.001 0.001];
  tspan1 = [-10 -1];
  iflag = 1;
  [t,z] = ode45(@(t,z) systode(t,z,iflag), tspan1, z0);
  T = vertcat(T,t);
  Z = vertcat(Z,z);
  tspan2 = [-1 0];
  iflag = 2;
  z0 = [z(end,1) z(end,2)];
  [t,z] = ode45(@(t,z) systode(t,z,iflag), tspan2, z0);
  T = vertcat(T,t);
  Z = vertcat(Z,z);
  tspan3 = [0 1];
  iflag = 3;
  z0 = [z(end,1) z(end,2)];
  [t,z] = ode45(@(t,z) systode(t,z,iflag), tspan3, z0);
  T = vertcat(T,t);
  Z = vertcat(Z,z);
  tspan4 = [1 10];
  iflag = 4;
  z0 = [z(end,1) z(end,2)];
  [t,z] = ode45(@(t,z) systode(t,z,iflag), tspan4, z0);
  T = vertcat(T,t);
  Z = vertcat(Z,z);
  figure
  plot(T,Z(:,1),'r');
end  

function dwdt = systode(t,w,iflag) 
  if iflag == 1
    f = 1/t;
    h = 0;
    beta=0.5+exp(t);
  elseif iflag == 2
    f = t*(3+2*t);
    h = 4*t^4+12*t^3+9*t^2+4*t+3;
    beta=0.5+exp(t);
  elseif iflag == 3
    f = t*(3-2*t);
    h = 4*t^4-12*t^3+9*t^2-4*t+3;
    beta=0.5+exp(-t);
  elseif iflag == 4
    f = 1/t;
    h = 0;
    beta = 0.5+exp(-t);
  end    
  
  dwdt=zeros(2,1);
  dwdt(1)=-f*w(1)+w(2);
  dwdt(2)=-beta*w(1)-f*w(2)+h*w(1)-f*w(1)^2;  
end 

[Cris19]

Thank you so very much, Torsten! Should I save each function 'main' and 'systode' separately? Or both in a single function file? I really don't know to handle these codes. And which are the commands for obtaining the plotting of the solution?

[Torsten]

Put the code in a file, name it main.m and load it into MATLAB. Run it and the first function will be plotted in the interval [-10:10]. The command is

figure
plot(T,Z(:,1),'r');

[Cris19]

Thank you a lot! It works very well.

Answer 2

[Cris19]

One more question, if possible...

I would also need to plot on the same graph the derivative of first component of the solution, dwdt(1). From a previous question posted on this forum, I learned that in the case when the coefficients are not piecewise defined, the code can be:

tspan = [-10 10];
z0=[0.001 0.001];
[t,z] = ode45(@(t,z) systode(t,z), tspan, z0);
for k = 1:numel(t)
    dwdt(:,k) = systode(t(k),z(k,:));
end
figure
plot(t,z(:,1),'b')
hold on
plot(t, dwdt(1,:), 'k')
hold off
legend('$x(t)$','$\dot{x}(t)$', 'Interpreter','latex', 'Location','best');

But in the present case, I can not handle it. Would you be kind to help me?

[Torsten]

function main
T = [];
Z = [];
dZ = [];
z0=[0.001 0.001];
tspan1 = [-10 -1];
iflag = 1;
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan1, z0);
dz = zeros(2,numel(t));
for i=1:numel(t)
    dz(:,i) = systode(t(i),z(i,:),iflag);
end  
T = vertcat(T,t);
Z = vertcat(Z,z);
dZ = vertcat(dZ,dz.');
tspan2 = [-1 0];
iflag = 2;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan2, z0);
dz = zeros(2,numel(t));
for i=1:numel(t)
    dz(:,i) = systode(t(i),z(i,:),iflag);
end  
T = vertcat(T,t);
Z = vertcat(Z,z);
dZ = vertcat(dZ,dz.');
tspan3 = [0 1];
iflag = 3;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan3, z0);
dz = zeros(2,numel(t));
for i=1:numel(t)
    dz(:,i) = systode(t(i),z(i,:),iflag);
end 
T = vertcat(T,t);
Z = vertcat(Z,z);
dZ = vertcat(dZ,dz.');
tspan4 = [1 10];
iflag = 4;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan4, z0);
dz = zeros(2,numel(t));
for i=1:numel(t)
    dz(:,i) = systode(t(i),z(i,:),iflag);
end 
T = vertcat(T,t);
Z = vertcat(Z,z);
dZ = vertcat(dZ,dz.');
figure
plot(T,Z(:,1),'r');
figure
plot(T,dZ(:,1),'g');
end  
function dwdt = systode(t,w,iflag) 
if iflag == 1
    f = 1/t;
    h = 0;
    beta=0.5+exp(t);
elseif iflag == 2
    f = t*(3+2*t);
    h = 4*t^4+12*t^3+9*t^2+4*t+3;
    beta=0.5+exp(t);
elseif iflag == 3
    f = t*(3-2*t);
    h = 4*t^4-12*t^3+9*t^2-4*t+3;
    beta=0.5+exp(-t);
elseif iflag == 4
    f = 1/t;
    h = 0;
    beta = 0.5+exp(-t);
end    

dwdt=zeros(2,1);
dwdt(1)=-f*w(1)+w(2);
dwdt(2)=-beta*w(1)-f*w(2)+h*w(1)-f*w(1)^2;  
end 

[Cris19]

Thank you so very much for the code, Torsten! I didn't know many things in Matlab and I'm starting to learn them.

[Torsten]

This is the more usual formulation for your ODE system.

Seems it was not absolutely necessary to interrupt integration at the points where the piecewise functions were glueed together.

function main
tspan = [-10 10];
z0 = [0.001 0.001];
[t,z] = ode15s(@(t,z) systode(t,z), tspan, z0);
for i=1:numel(t)
    dz(:,i) = systode(t(i),z(i,:));     
end
figure
plot(t,z(:,1),'r')
figure
plot(t,dz(1,:).','g')
end   

function dwdt = systode(t,w)     
if abs(t) < 1 
    f = t*(3-2*abs(t));
else
    f = 1/t;
end

if t>=0 && t < 1 
    h=4*t^4-12*t^3+9*t^2-4*t+3;
elseif t>-1 && t <= 0
    h=4*t^4+12*t^3+9*t^2+4*t+3;
else 
    h=0;
end  

beta=0.5+exp(-abs(t));
dwdt=zeros(2,1);
dwdt(1)=-f*w(1)+w(2);
dwdt(2)=-beta*w(1)-f*w(2)+h*w(1)-f*w(1)^2;
end    

[Cris19]

Yes, indeed, it seems to be easier this way. I think it is possible, since the functions f, h, β are continuous (and also differentiable) on the whole [-10,10].