# How to solve a 2nd-order ODE system with space-dependent variable with ODE45?

1 view (last 30 days)
Ying Wu on 9 Aug 2021
Commented: Ying Wu on 13 Aug 2021
Hi, I use ode45 function to solve the following 2nd-order ODE system:  where X(t) and Y(t) are function of time, and the parameter P is dependent on the first dirivative of Y(t).
I have 15 scattered points for P and dY:
dY = [0, 0.0026, 0.0082, 0.0138, 0.0194, 0.0250, 0.0306, 0.0362, 0.0418, 0.0475, 0.0532, 0.0589, 0.0647, 0.0705, 0.0763]
P = [0, 0.0423, 0.1534, 0.2336, 0.2915, 0.3628, 0.4166, 0.4587, 0.5572, 0.8002, 0.8204, 0.5559, 0.3184, 0.1405, -0.0603]
My main steps include:
1. I use cubic spline interpolation to fit these points and obtain a piecewise function for P(dY): pp=csapi(dY,P).
2. Like the time-dependent variable in Matlab ode45 help, I define a vector of dY and obtain the corresponding P(dY).
3. Treat (dY,P) as inpus of odefun and use interp1 to update P with dY at each time step during the numerical computation.
The key code are shown below:
% cubic spline interpolation
csi = csapi(dY,P);
dY = [0:0.0001:0.1]';
P = fnval(csi,dy);
% ode45
[time,sol] = ode45(@(t,Y) odefun(t,Y,dY,P),[0 1000],[0,0,0.01,0);
function dYdt = odefun(t,Y,dY,P)
CCFy = interp1(ddY,P,abs(Y(4)));
% Y(1) = X
% Y(2) = X'
% Y(3) = Y
% Y(4) = Y'
dYdt = zeros(4,1);
dYdt(1) = Y(2);
dYdt(2) = 0.51*(0.03*Y(1) + 0.0025*P - 0.045*Y(4) - Y(3)) + 0.887*Y(4) + ...
1.1642*(1-1175.51*Y(1)^2)*Y(2)- 1.774*Y(1);
dYdt(3) = Y(4);
dYdt(4) = 0.03*Y(1) + 0.0025*P - 0.045*Y(4) - Y(3);
end
But the results are incorrect. Could anyone please give me some help? Thanks!

Alan Stevens on 9 Aug 2021
More like this perhaps
dY = [0, 0.0026, 0.0082, 0.0138, 0.0194, 0.0250, 0.0306, 0.0362, 0.0418, 0.0475, 0.0532, 0.0589, 0.0647, 0.0705, 0.0763];
P = [0, 0.0423, 0.1534, 0.2336, 0.2915, 0.3628, 0.4166, 0.4587, 0.5572, 0.8002, 0.8204, 0.5559, 0.3184, 0.1405, -0.0603];
tspan = [0 1000];
ic = [0,0.01,0,0];
[time,sol] = ode45(@(t,Y) odefun(t,Y,dY,P),tspan,ic);
subplot(2,1,1)
plot(time,sol(:,1)),grid
xlabel('time'),ylabel('Y')
subplot(2,1,2)
plot(time,sol(:,3)),grid
xlabel('time'),ylabel('X') function dYdt = odefun(~,Y,dY,P)
% Y(1) = Y
% Y(2) = Y'
% Y(3) = X
% Y(4) = X'
p = @(dydt) spline(dY,P,dydt);
dYdt = zeros(4,1);
dYdt(1) = Y(2);
dYdt(2) = 0.03*Y(3) + 0.0025*p(Y(2)) -0.045*Y(2);
dYdt(3) = Y(4);
dYdt(4) = 0.51*dYdt(2) + 0.887*Y(2) - 1.774*Y(3) + ...
1.1642*(1-1175.51*Y(3)^2)*Y(4);
end
Ying Wu on 13 Aug 2021
Thanks again for your suggestion!

R2020b

### Community Treasure Hunt

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

Start Hunting!