# How to use ODE45 for a coupled system of differential equations at a specific time point

3 views (last 30 days)

Show older comments

Hey everyone. I have been using ode45 for a while and it was smooth for me. However, I am having trouble doing it for this system because its a bit complicated. I had multiple attempts fail before giving up. I hope someone can help.

I have the system shown below in the picture, its a system of three differential equations. Wx, Wy, and Wz are constants that I define, and the initial conditions for theta, phi, and psi, are also constants that I define. I want to find the solution for theta, phi and psi at a specific point of time only. Example, time = 2.2 seconds. Note that the dots above the variables indicate the derivative, so x_dot is dx/dt.

##### 3 Comments

Sam Chak
on 6 Apr 2023

### Accepted Answer

Jon
on 6 Apr 2023

So, expanding on @Davide Masiello comment, if all you want is the final solution at t = 2.2s you could wrap your call to ode45 using something like this (numerical values are just for illustrative purposes)

theta0 = pi/10;

phi0= pi/3;

psi0 = pi/9;

w = randn(3,1);

tFinal = 2.2

[theta,phi,psi] = solveSys(tFinal,theta0,phi0,psi0,w)

where solveSys is a function that you define (store as its own .m file )

function [theta,phi,psi] = solveSys(tFinal,theta0,phi0,psi0,w)

% solve system of odes to get final value of angles theta,phi and psi for

% specified w = [wx,wy,wz]

tspan = [0,tFinal];

x0 = [theta0,phi0,psi0]; % initial conditions

% define some local variables for readability

wx = w(1);

wy = w(2);

wz = w(3);

[~,x] = ode45(@fun,tspan,x0); % just get back angle values, not interested in times

% just return the final values

theta = x(end,1);

phi = x(end,2);

psi = x(end,3);

function xdot = fun(~,x)

% calculates state derivatives, time is not needed, but needs to be

% included in signature, so use ~

% use nested function so parameters wx,wy,wz will be

% available

% define states for clearer readability

theta = x(1);

phi = x(2);

psi = x(3);

phiDot= wx + tan(theta)*sin(phi)*wy + tan(theta)*cos(phi)*wz;

thetaDot = cos(phi)*wy - sin(phi)*wz;

psiDot = sec(theta)*sin(phi)*wy + sec(theta)*cos(phi)*wz;

xdot = [phiDot;thetaDot;psiDot];

end

end

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!