It works for me (I've commented out the fprint statements in the code below - but it works with them in!):
global CDin psi theta1 theta2 a b
% fileID = fopen('Project_ode45.txt','w');
% Define Parameters for the problem
CDin = 4; %mol/L of Co2 inlet
psi = 0.1809; %psi=v/(V*k2)
theta1 = 0.0879*CDin^2; %(k1/k2)*Cdin^2 with
theta2 = (3.587/10000)*CDin; %k3/k2
a = 3; %N/C need to specify
b = 0.6; %H/C need to specify
y0 = [a;0;0;1;0;b]; %only co2 is 1, we can change a and b
tspan = [0 2];
[T,Y] = ode45(@func,tspan,y0);
dummy = [T,Y]; % Store output into a dummy matrix28
plot(dummy(:,1),dummy(:,2),dummy(:,1),dummy(:,3),dummy(:,1),dummy(:,4),dummy(:,1 ),dummy(:,5),dummy(:,1),dummy(:,6),dummy(:,1),dummy(:,7)) % Plot the solution
% fprintf( ' Time y1 y2 y3 y4 y5 y6\n'); % Print output headings
% fprintf(fileID,' Time y1 y2 y3 y4 y5 y6\r\n'); % Print output headings
% fprintf( '%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n',dummy');
% % Print to screen
% fprintf(fileID,'%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\r\n',dummy');
% Print to file
function f = func(t,y)
% function func to input the RHS of IVP's
global psi theta1 theta2 a b
f = zeros(length(y),1)+t*0; % Reset RHS functions
f(1) = psi*(a-y(1))-2*theta1*y(1)^2*y(4)+theta2*y(5)^2;
f(2) = theta2*y(5)^2-psi*y(2);
f(3) = theta1*y(1)^2*y(4)-y(3)-psi*y(3);
f(4) = psi*(1-y(4))-theta1*y(1)^2*y(4);
f(5) = y(3)-2*theta2*y(5)^2-psi*y(5);
f(6) = psi*(b-y(6))+y(3);
end
