2nd order differential eqn for Windkessel model

I am struggling to solve this 2nd order ODE in matlab for a 4WK Parallel Windkessel model, this is what I have so far for my input and outputs for a 2WK model, a 3WK model (attached).

I am trying to solve a 2nd order ODE:

I''(R*L*C*r) + I'(L*(R+r)) +I*(R*r) = Y"(L*C*R)+Y'(C*R*r +L)+ Y*r

I is :

I=@(t)I0*sin((pi*t)/Ts).^2.*(t<=Ts); %input current flow

t=0:0.001:Tc;

Initial conditions for Y(0)= 80, the rest of the variables are known constants.

Thanks in advance!

[Torsten]

You have a second-order ODE in Y. So you need to solve a system of equations for Y and Y'.

[Nang Su Lwin]

Yes, I think the errors are stemming from my function handles of I though, since most examples I have seen just have constants rather than functions of time. The complications are coming from the 2nd order ODE, and the differentials in I that I am trying to solve. I was able to solve for Y fine with first order ODE's as seen in my code.

[Torsten]

I was able to solve for Y fine with first order ODE's as seen in my code.

Then you must have used an ODE different from the one you posted:

I''(R*L*C*r) + I'(L*(R+r)) +I*(R*r) = Y"(L*C*R)+Y'(C*R*r +L)+ Y*r

[Nang Su Lwin]

Yes, there are three different ODE's being solved in the code. I was able to solve for the first two ODE's because they were only first order ODE's. For, this last ODE that I posted I can't find a way to have a system of equations to solve for it since I have functions in time for I and the function is a piecewise function (see initial post where it is a sin^2(t) only for <Ts and zero for the rest of the time period). Ts, and the time period (0.3333, 0.83333) are called in the code

[Torsten]

You have to solve

y1' = y2
y2' = (I''*(R*L*C*r) + I'*(L*(R+r)) +I*(R*r) - (y2*(C*R*r +L)+ y1*r))/(L*C*R)

with

y1(0) = 80 and y2(0) = ? 

(y1 = Y, y2 = Y')

[Nang Su Lwin]

I think y2(0) = 0

How would I call these system of equations and include the I functions of time?

[Nang Su Lwin]

I tried to simplify the equations as much as possible to help me not go crazy as I tried to solve this equation (I substituted in all the constants so I didn't have to deal with so many variables). The problem is that I don't know how to implement the fact that my Input function I,I', and I'' (now to the right of the equation) doesn't consider that they should only output when t<0.33333:

syms y(t)

Dy = diff(y,t);

D2u = diff(y,t,2);

ode= 0.05*y+0.055*diff(y,t,1)+ +0.005*diff(y,t,2) == 15.75*pi*cos(3*pi*t)*sin(3*pi*t)+22.5*pi^2*cos(3*pi*t)^2 - 22.5*pi^2*sin(3*pi*t)^2 +25*sin(3*pi*t).^2;

cond1 = y(0) == 80;

cond2 = Dy(0) == 0;

conds = [cond1 cond2];

ySol(t) = dsolve(ode, conds);

ySol= simplify(ySol)

t=0:0.001:60/72;

plot(t,ySol(t),'g')

So I get something like this:

Instead of what I'm expecting:

Answer 1

[Torsten]

I = @(t) I0*sin((pi*t)/Ts).^2.*(t<=Ts); %input current flow
Idot = @(t) I0*2*sin(pi*t/Ts).*cos(pi*t/Ts)*pi/Ts.*(t<=Ts);
Idotdot = @(t) I0*2*(cos(pi*t/Ts).^2 - sin(pi*t/Ts).^2)*(pi/Ts)^2.*(t<=Ts);
fun = @(t,y)[y(2);(Idotdot(t)*(R*L*C*r)+Idot(t)*(L*(R+r))+I(t)*(R*r) - (y(2)*(C*R*r+L)+y(1)*r))/(L*C*R)];
y0 = [80 ;0];
tspan = [0 5/6];
[T,Y] = ode45(fun,tspan,y0)

[Nang Su Lwin]

Thank you! That was a lot simpler than I expected, issue was not knowing the "y(2);" syntax!

Much appreciated!!

[Nang Su Lwin]

Just for others who might read this post later and are trying to also plot some Windkessel responses, this is the 4Windkessel code:

clear all
clc   
I0= 500; % maximum flow
Tc=60/72;   % heart period
Ts=(2/5*Tc);   % time in systole
P_ss=80; % diastolic pressure
R= 1;
C=1;
R1=0.05;    
L=0.005;
I=@(t)I0*sin((pi*t)/Ts).^2.*(t<=Ts); %input current flow
Idot = @(t)I0*2*sin(pi*t/Ts).*cos(pi*t/Ts)*pi/Ts.*(t<=Ts);
Idotdot = @(t) I0*2*(cos(pi*t/Ts).^2 - sin(pi*t/Ts).^2)*(pi/Ts)^2.*(t<=Ts);
fun = @(t,y)[y(2);(Idotdot(t)*(R*L*C*R1)+Idot(t)*(L*(R+R1))+I(t)*(R*R1) - (y(2)*(C*R*R1+L)+y(1)*R1))/(L*C*R)];
y0 = [80 ;0];
tspan = [0 60/72];
[T,Y] = ode45(fun,tspan,y0);
plot(T,Y(:,1),'g')