How to graph a complex system of ODEs?

Any suggestions on how to graph this system?

(1) dU/dt = C1 * ( 1 - C2*U)

(2) dV/dt = C2 * dU/dt + C3 * (V - W)

(3) dW/dt = C4 * (V - W)

Where the Cs are all unique constants

I tried using ode45 to create arrays but am having trouble figuring out how to weave together several differential MATLAB functions that all depend on one another especiallly given how equation 2 depends on equation 1 and equation 2 and 3 are intertwined with shared variables.

This is a system of equations describing an isothermal kinetics PFR reactor w/ Heat Exchanger in terms of simplified mathematical variables (rather than the many constants and X, T, Ta, and V variables). Normally PolyMath is used but that program is very buggy.

Answer 1

[James Tursa]

You have three 1st order equations, so you need a 3-element state vector for this. To keep the same nomenclature as the MATLAB docs, I will use y for this state vector. The definitions would be

y(1) = U

y(2) = V

y(3) = W

Then given a 3-element state vector y, write a derivative function for ydot. E.g.,

function ydot = myderivative(t,y,C1,C2,C3,C4)
ydot = zeros(size(y));
ydot(1) =
ydot(2) =
ydot(3) =
end

You write code for the three lines above. Wherever you see a U you would use y(1). Wherever you see a V you would use y(2). Wherever you would see a W you would use y(3). And dU/dt is simply ydot(1), be sure to calculate that first so you can use it in ydot(2). So write that code. Then your ode45( ) call would look something like:

[T,Y] = ode45(@(t,y)myderivative(t,y,C1,C2,C3,C4,tspan,y0);

[David Goodmanson]

Just a quick comment on James' answer. By default, building up the elements of ydot one by one gives a row vector. But this will error out because the output of the ydot function has to be a column vector. So within the function you can either preallocate ydot as zeros(3,1) first or use ydot = ydot' afterwards. ( Or build it up as ydot = [expression1;expression2;expression3] which can lead to a pretty dense line of code).

[James Tursa]

Good comment, but note that ydot = zeros(size(y)) accomplishes this since the incoming y is a column vector.

[Steven Lord]

To make the code in myderivative look even more similar to the equations, you can define some intermediate variables.

function ydot = myderivative(t,y,C1,C2,C3,C4)
U = y(1);
V = y(2);
W = y(3);
dUdt = % write it terms of t, U, V, W, etc.
dVdt =
dWdt = 
ydot = [dUdt; dVdt; dWdt];
end