How to graph two solutions in one continuous plot?

I am trying to model the solution to a system of ODEs (written in another file). This file runs and produces a plot, but I can't seem to join the plot for the two different solutions. (There are two ode45 solvers here.) I want the first one to go from t=0 to t=5, the second from t=5 to t=120.
I have tried defining piecewise functions, tried defining two different lin spaces, and it's not working.
Any suggestions would be helpful!
Thank you!
Below is my code:
clear all
clc
axis manual;
close all;
global lamdax Dx beta q Dy eta k u lamdaz c C1 K1 Dw p C2 K2 rho Dz f Cm C3 K3 Dm
lamdax=100;
lamdaz = 600/24;
Dx=.1;
Dy = .5;
Dw = 1/24;
Dz = 1/40;
Dm = 1/1500;
u = 3;
beta = .00061;
eta = .01;
k = 100;
q = 1/30;
c = .006;
Cm = 1.6;
f = 4;
p = .6;
rho = .06;
K1 = 10;
K2 = 10;
K3 = 10;
C1 = .1;
C2 = .1;
C3 = .1;
%parameter values
%solving the ODE @Bliss_ode for time span [0, 1500] and initial
%conditions x1(0)=x2(0)=x3(0)=0.5
[t,y] = ode45(@measles_ode, [0 150], [lamdax/Dx; .001; 300; lamdaz/Dw; .001; .001]);
%[t,y] = ode45(@measles_ode, [5, 150], [993.5; 2.75; 300; 595.7; 7.8; .57]);
%, options);%,[],parfit);
%%plot of x1 versus time
subplot(1,6,1)
plot(t, y(:,1), 'r', 'LineWidth',2);
%%plot of x2 versus time
subplot(1,6,2)
plot(t, y(:,2), 'b', 'LineWidth',2);
%plot of x3 versus time
subplot(1,6,3)
plot(t, y(:,3), 'color',[0.9100 0.4100 0.1700], 'LineWidth',2);
%xlabel('time in days')
%%ylabel('Virus')
subplot(1,6,4)
plot(t, y(:,4), 'k', 'LineWidth',2);
subplot(1,6,5)
plot(t, y(:,5), 'y', 'LineWidth',2);
subplot(1,6,6)
plot(t, y(:,6), 'g', 'LineWidth',2);

 Accepted Answer

EDITED
clear all
clc
axis manual;
close all;
global lamdax Dx beta q Dy eta k u lamdaz c C1 K1 Dw p C2 K2 rho Dz f Cm C3 K3 Dm
lamdax=100;
lamdaz = 600/24;
Dx=.1;
Dy = .5;
Dw = 1/24;
Dz = 1/40;
Dm = 1/1500;
u = 3;
beta = .00061;
eta = .01;
k = 100;
q = 1/30;
c = .006;
Cm = 1.6;
f = 4;
p = .6;
rho = .06;
K1 = 10;
K2 = 10;
K3 = 10;
C1 = .1;
C2 = .1;
C3 = .1;
%parameter values
%solving the ODE @Bliss_ode for time span [0, 1500] and initial
%conditions x1(0)=x2(0)=x3(0)=0.5
[t,y] = ode45(@measles_ode, [0 150], [lamdax/Dx; .001; 300; lamdaz/Dw; .001; .001]);
[t1,y1] = ode45(@measles_ode, [5, 150], [993.5; 2.75; 300; 595.7; 7.8; .57]);
%, options);%,[],parfit);
%%plot of x1 versus time
subplot(6,1,1)
plot(t, y(:,1), 'r', 'LineWidth',2);
hold on
plot(t1, y1(:,1), 'b', 'LineWidth',2);
%%plot of x2 versus time
subplot(6,1,2)
plot(t, y(:,2), 'b', 'LineWidth',2);
hold on
plot(t1, y1(:,2), 'r', 'LineWidth',2);
%plot of x3 versus time
subplot(6,1,3)
plot(t, y(:,3), 'color',[0.9100 0.4100 0.1700], 'LineWidth',2);
hold on
plot(t1, y1(:,3), 'b', 'LineWidth',2);
%xlabel('time in days')
%%ylabel('Virus')
subplot(6,1,4)
plot(t, y(:,4), 'k', 'LineWidth',2);
hold on
plot(t1, y1(:,4), 'y', 'LineWidth',2);
subplot(6,1,5)
plot(t, y(:,5), 'y', 'LineWidth',2);
hold on
plot(t1, y1(:,5), 'g', 'LineWidth',2);
subplot(6,1,6)
plot(t, y(:,6), 'g', 'LineWidth',2);
hold on
plot(t1, y1(:,6), 'b', 'LineWidth',2);
function dx = measles_ode(t,x)
global lamdax Dx beta q Dy eta k u lamdaz c C1 K1 Dw p C2 K2 rho Dz f Cm C3 K3 Dm
dx = zeros(6,1);
dx(1)= lamdax-Dx*x(1)-beta*q*x(1)*x(3);
dx(2)= beta*q*x(1)*x(3)-Dy*x(2)-eta*x(2)*x(5);
dx(3)= k*x(2)-u*x(3)-beta*q*x(3)*x(1);
dx(4)= lamdaz-((c*q*x(4)*x(3))/(C1*q*x(3)+K1))-Dw*x(4);
dx(5)= ((c*q*x(4)*x(3))/(C1*q*x(3)+K1))+((p*q*x(3)*x(5))/(C2*q*x(3)+K2))-(rho+Dz)*x(5)...
+((f*Cm*q*x(3)*x(6))/(C3*q*x(3)+K3));
dx(6)= (rho*x(5))-(Dm*x(6))-((Cm*q*x(3)*x(6))/(C3*q*x(3)+K3));
end
Screen Shot 2018-11-23 at 9.13.18 AM.png

6 Comments

Hello! Thanks for the input.
I've tried looking at examples similar to this. I'm not sure how to assign y1 and y2.
I've tried:
y1 = [t,y] = ode45(@measles_ode, [0 150], [lamdax/Dx; .001; 300; lamdaz/Dw; .001; .001]);
y2 = [t,y] = ode45(@measles_ode, [5, 150], [993.5; 2.75; 300; 595.7; 7.8; .57]);
which gives me an error saying that I can't use =, I need ==.
so I tried y1 = [t,y] ==
as well as y1== [t,y] =
neither of which worked. I've also tried replacing [t,y] with just y1 or y2, as well as [t, y1] and [t, y2] which also didn't work
I've tried all of these things, but none of them worked. I'm relatively new to MatLab, so I apologize. If you could show me an example of how I should define y1 and y2, I would appreciate it!
Thank you!
Upload your function and explain which solutions you want to stitch
Ideally, I would like to stitch together these two solutions:
(1) [t,y] = ode45(@measles_ode, [0 150], [lamdax/Dx; .001; 300; lamdaz/Dw; .001; .001]);
(2) [t,y] = ode45(@measles_ode, [5, 150], [993.5; 2.75; 300; 595.7; 7.8; .57]);
I'd like (1) to be on the interval from [0, 5] and (2) to be on the interval from [5, 120].
The thought process behind it being that the first one models the measles virus vaccination in the first 5 years, and then there's a booster shot at year 5, which is modeled by the second one. (Uses the values at the end of year 5 for the initial conditions in the second solver, except the therapy which stays the same - a 'booster' of 300.)
Like I said, I'm very inexperienced with MATLAB so this is just an intuition on how I might show this scenario. I know the first one is correct, I'm just not sure about showing the 'booster' portion, and how to fit the two together on the same graph.
The ODE system, if it helps is the following:
function dx = measles_ode(t,x)
global lamdax Dx beta q Dy eta k u lamdaz c C1 K1 Dw p C2 K2 rho Dz f Cm C3 K3 Dm
dx = zeros(6,1);
dx(1)= lamdax-Dx*x(1)-beta*q*x(1)*x(3);
dx(2)= beta*q*x(1)*x(3)-Dy*x(2)-eta*x(2)*x(5);
dx(3)= k*x(2)-u*x(3)-beta*q*x(3)*x(1);
dx(4)= lamdaz-((c*q*x(4)*x(3))/(C1*q*x(3)+K1))-Dw*x(4);
dx(5)= ((c*q*x(4)*x(3))/(C1*q*x(3)+K1))+((p*q*x(3)*x(5))/(C2*q*x(3)+K2))-(rho+Dz)*x(5)...
+((f*Cm*q*x(3)*x(6))/(C3*q*x(3)+K3));
dx(6)= (rho*x(5))-(Dm*x(6))-((Cm*q*x(3)*x(6))/(C3*q*x(3)+K3));
Thank you so much madhan ravi! You're the best! That worked!!
Anytime :) , make sure to accept the answer if it helped

Sign in to comment.

More Answers (0)

Categories

Find more on Mathematics in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!