The error you made was in not putting them all in the same system. You did code them correctly otherwise, so all I needed to do was to copy your code to the ‘cprime’ matrix I created, subscript the concentrations, and provide the correct initial conditions vector. You have to plot them using the subplot function or you won’t be able to see the changes.
This works:
a = 1;
b = (2E-3);
d = (4E-4);
f = (4E-8);
C10 = .063;
C20 = 700;
C30 = 5.6;
cprime = @(t,C) [(-a*(C(1).^3)) + (b*(C(3))) + (d*(C(3).^4)) - f*(C(2).^4).*(C(1).^4); (d*(C(3).^4)) - f*(C(2).^4).*(C(1).^4); (a*(C(1).^3)) - (b*C(3)) - (d*(C(3).^4)) + f*(C(2).^4).*(C(1).^4)]
options = odeset('RelTol',10^-7,'AbsTol',10^-9);
[X,C] = ode45(cprime, [0 1], [C10; C20; C30], options);
figure(1)
subplot(3,1,1)
plot(X, C(:,1))
ylabel('\itC_1\rm')
grid
subplot(3,1,2)
plot(X, C(:,2))
ylabel('\itC_2\rm')
grid
subplot(3,1,3)
plot(X, C(:,3))
grid
xlabel('\itx\rm')
ylabel('\itC_3\rm')
I experimented with ode15s as well. It gives the same results as ode45.
