Creating a 3d plot of the solutions to an ODE as the time a switch condition activates varies

Hi!

I am using a function with a subfunction to produce a graph of a piecewise smooth variable in an ODE. Smoothness is lost when t=t_1. Code for this is shown at the bottom of the question.

I'd like code to produce a 3d graph showing how different solutions look as t_1 varies. This would have the three axes t, t_1 and xa(:,1). Currently what I am stuck on is I am not sure how I could include t_1 as a a variable that is fed into ydot without creating confusing when ode45 calls it.

function L2_ode45
%a function to call ode45 and draw the graph
  [t,xa] = ode45(@fu,[0 500],[0 0]);
  plot(t,xa(:,1))
  grid on
end

function ydot = fu(t,x)
%the function passed to ode45
%give values to parameters
t_1  = 200
k_a1 = 10^(-4);
k_a2 = 2*10^(-2);
C    = 1;
R    = 10;
k_d1 = 10^(-3);
k_d2 = 10^(-3);
%set up the switch
  if t < t_1
    z = C*(R-x(1)-x(2));
  else
    z = 0;
  end
%set up the ROC function that will be output to ode45
ydot = [  k_a1*z-k_d1*x(1);...
          k_a2*z-k_d2*x(2);];
end

Answer 1

[Star Strider]

I am not certain what you want to do.

Try this:

function L2_ode45
    tspan = linspace(0, 500, 50);
    t1 = 100:100:500;
    for k = 1:numel(t1)
        [t{k},xa{k}] = ode45(@(t,x)fu(t,x,t1(k)),tspan,[0 0]);
    end

q1 = size(t)
tp = cell2mat(t);
xap = cell2mat(xa);
xapm = xap(:,1:2:end);

figure
mesh(t1, tspan, xapm)
grid on
xlabel('t_1')
ylabel('t')
end

function ydot = fu(t,x,t_1)
%give values to parameters
% t_1  = 200
k_a1 = 10^(-4);
k_a2 = 2*10^(-2);
C    = 1;
R    = 10;
k_d1 = 10^(-3);
k_d2 = 10^(-3);
%set up the switch
  if t < t_1
    z = C*(R-x(1)-x(2));
  else
    z = 0;
  end
%set up the ROC function that will be output to ode45
ydot = [  k_a1*z-k_d1*x(1);...
          k_a2*z-k_d2*x(2);];
end

If you want to interrupt the integration to substitute some new value, see the documentation on ODE Event Location.

EDIT —

I forgot to post the plot!

Here it is —