If you have the right hand side of your ODE as an n-by-1 vector, then you can just assign the output of the ode function to this vector.
As an example, below I have implemented the simulation of a rotating rigid body with a constant angular velocity about the body x-axis. Note that the equations of motion of a rotating rigid body, as expressed in the body frame, are
R'(t) = R(t)*Omega(t);
I*omega'(t) = I*omega(t) x omega(t)
where R is the rotation matrix, Omega is the 3-by-3 skew-symmetric angular velocity matrix, omega is the representation of the angular velocity as a 3-vector (a vector in R^3), and I is the moment of inertia of the body expressed in the body frame. The primes denote time differentiation and "x" denotes the cross product in R^3.
In particular, the function "odefun" is doing what you are asking. It creates the vector field "f" and assigns it as the output of "odefun".
function odeSystemSample()
% Simulation of a rotating rigid-body
%%Set-up the problem and solve
par.I = diag([1,2,3]); % Inertia of the body
q0 = [reshape(eye(3),9,[]); 1; 0; 0]; % Initial condition, body is
% rolling with constant angular
% velocity
ti = 0; tf = 2;
opt = odeset('AbsTol',1.0e-07,'RelTol',1.0e-07);
[t,q] = ode45( @odefun, [ti:0.001:tf], q0 ,opt, par);
end
function dq = odefun(t,q,par)
% This is how you assign an n-by-1 vector to the output of the function
% that defines the differential equations.
R = reshape(q(1:9),3,3);
w = q(10:12);
dq = zeros(12,1);
dq(1:9) = reshape(R*phi(w,0),9,[]); % This creates a 9-by-1 vector
dq(10:12) = par.I \ cross(par.I*w, w); % This is a 3-by-1 vector
end
function x = phi(y,flag)
% This function is problem specific; you do not need to pay attention to it.
% This is the isomorphism between the Lie algebra so(3) and R^3.
% If the flag is on, then the function goes from so(3) to R^3.
% If the flag is off, then the function goes from R^3 to so(3).
if flag
x = [-y(2,3); y(1,3); -y(1,2)];
else
x = [0, -y(3), y(2); y(3), 0, -y(1); -y(2), y(1), 0];
end
end
