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

9 views (last 30 days)
Harry on 15 Jan 2020
Edited: Star Strider on 15 Jan 2020
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

Star Strider on 15 Jan 2020
Edited: Star Strider on 15 Jan 2020
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 — 