I have solved a second order differential equation $\ddot R=f(t,R,\dot R)$ with a piecewise function $$P(t)=\begin{cases}P_0 & \text{ if } 0<t\leq40\times10^{-6} \\ \frac{(P_0+P_{min})}{2}+\frac{(P_0-P_{min})}{2}\cos(2\pi ft) & \text{ if } 40\times10^{-6}<t\leq60\times10^{-6} \\ P_{min} & \text{ if } 60\times10^{-6}<t\leq100\times10^{-6} \\ \frac{(P_0+P_{min})}{2}+\frac{(P_0-P_{min})}{2}\cos(2\pi ft) & \text{ if } 100\times10^{-6}<t\leq120\times10^{-6} \\ P_0 & \text{ if } t>120\times10^{-6} \end{cases}$$I would like to plot, on the same graph with two $y$ axes (I guess with the plotyy command) both $t$ versus $R$ and $t$ versus $P(t)$. At the moment I am only able to plot $t$ versus $R$ plot(t,R(:,1)).function overshoot()R0 = 15e-6;tspan = [0 160e-6];[t,R] = ode45(@(t,R) DE(t,R,R0), tspan, [R0,0]);t = t * 1e6;R = R * 1e6;plot(t,R(:,1))end%function Rdot = DE(t,R,R0)S = 0.073;rho = 998;mi = 1.005e-3; % kg / (m * s)P0 = 101325;Pvap = 2329.6;Pmin = 1800;f = 1e6 / 40;%Rdot = zeros(2,1);Rdot(1) = R(2);Rdot(2) = -1.5 * R(2) * R(2) / R(1) + 1 / (R(1) * rho) *...    (Pvap - P(t,f,P0,Pmin) + (P0 - Pvap + 2 * S / R0) *...    (R0 / R(1))^3 - 2 * S / R(1) - 4 * mi * R(2) / R(1));end%function fval = P(t,f,P0,Pmin)if (t <= 40e-6)    fval = 101325;elseif (t > 40e-6) && (t <= 60e-6)    fval = (P0 + Pmin) / 2 + (P0 - Pmin) / 2 * cos(2 * pi * f * t);elseif (t > 60e-6) && (t <= 100e-6)    fval = 1800;elseif (t > 100e-6) && (t <= 120e-6)    fval = (P0 + Pmin) / 2 + (P0 - Pmin) / 2 * cos(2 * pi * f * t);else    fval = 101325;endend
ODE in Matlab: how to plot result and piecewise function
matlab
This can achieved by making the parameters P0 Pmin f global.function overshoot()global P0 Pmin fR0 = 15e-6;tspan = [0 160e-6][t,R] = ode45(@(t,R) DE(t,R,R0), tspan, [R0,0]);%t = t * 1e6;%R = R * 1e6;figure(1)plot(t,R(:,1),'r--')hold onplot(t,P(t,f,P0,Pmin),'o');end  function Rdot = DE(t,R,R0)global P0 Pmin fS = 0.073;rho = 998;mi = 1.005e-3; % kg / (m * s)P0 = 101325;Pvap = 2329.6;Pmin = 1800;f = 1e6 / 40;%Rdot = zeros(2,1);Rdot(1) = R(2);Rdot(2) = -1.5 * R(2) * R(2) / R(1) + 1 / (R(1) * rho) *...    (Pvap - P(t,f,P0,Pmin) + (P0 - Pvap + 2 * S / R0) *...    (R0 / R(1))^3 - 2 * S / R(1) - 4 * mi * R(2) / R(1));endfunction fval = P(t,f,P0,Pmin)for i=1:length(t)if (t(i) <= 40e-6)    fval(i) = 101325;elseif (t(i) > 40e-6) && (t(i) <= 60e-6)    fval(i) = (P0 + Pmin) / 2 + (P0 - Pmin) / 2 * cos(2 * pi * f * t(i));elseif (t(i) > 60e-6) && (t(i) <= 100e-6)    fval(i) = 1800;elseif (t(i) > 100e-6) && (t(i) <= 120e-6)    fval(i) = (P0 + Pmin) / 2 + (P0 - Pmin) / 2 * cos(2 * pi * f * t(i));else    fval(i) = 101325;endend    endThe output plot which now has two $y$ axes, i.e, $R(t)$ and $P(t)$. But the plot doesn't make any sense because the maximum for $R(t)$ is approximately  $8\times10^{-5}$ and for $P(t)$ is $\approx 10\times10^{4}$.Anyways the desired combined plot looks like thisI hope this helps.