Function error in ODEs RK4 method, solving 6 unknowns

If I run this code, I get the following error.

Unrecognized function or variable 'func'.

Error in odearguments (line 92)

f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 104)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

% Solves a system of IVP's using ode45

% Term Project (Dimensionless)

clear, clc

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

The following equations are need to be solved. y1 to y6 are the variables that are need to be find.

Answer 1

[Alan Stevens]

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

[Common]

Thanks for your quick response

I still get the same error

I afraid that whether It only works with licenced version. I'm using Matlab free online basic version.

[VBBV]

Don't run the code in command window. Copy the code to a script file and save it using a valid filename, then run it from editor window using the Run button (green color)

[Common]

Thanks, It works. Im a beginner, pls dont mind