ODE45 - seperate function from commandfile

I'm calculating some chemical reactions using the code below but I would like to bring out the "incode" function concentration to a seperate functionfile. I've made a try in the code at the bottom.

s = 1;
q = 1; 
w = 0.1610;

% define y = [α(t) β(t) γ(t)] = [y(1) y(2) y(3)], then

concentration = @(t,y) [...
    s * (y(2) - y(1)*(y(2)+1) - q*y(1)^2)   % α'(t)
    (-y(2) - y(1)*y(2) + y(3))/s            % β'(t)
    w * (y(1)-y(3))                         % γ'(t)
];

% initial value
y0 = [30 1 30];

% time span 
t_span = [0 10];

% solve numerically
[t, y] = ode45(concentration, t_span, y0);

plot(t,y);

Function draft

function concentration = func(t,y) 

concentration=[
    s * (y(2) - y(1)*(y(2)+1) - q*y(1)^2)   % α'(t)
    (-y(2) - y(1)*y(2) + y(3))/s            % β'(t)
    w * (y(1)-y(3))                         % γ'(t)
];

I'm a bit lost and would appreciate some help to do this.

[Jan]

Uh, the lovely "urgent" tag. See http://www.mathworks.com/matlabcentral/answers/29922-why-your-question-is-not-urgent-or-an-emergency.

Do you get an error message, when you run this? Or what else is the problem?

[Matt Kindig]

Well for one, ode45 takes a function as its first argument, rather than a variable. As your function is currently defined, the function is named 'func', without an output of 'concentration'. You probably want the function itself to be called concentration, such as:

function c = concentration(t,y)

c=[
    s * (y(2) - y(1)*(y(2)+1) - q*y(1)^2)   % α'(t)
    (-y(2) - y(1)*y(2) + y(3))/s            % β'(t)
    w * (y(1)-y(3))                         % γ'(t)
];

[Jan]

@fxo: Not funny. People have died.

[Jan]

@fxo: Please note, that the "why is your question not urgent" thread is the one with the 2nd most votes in this forum. See also: http://www.mathworks.com/matlabcentral/answers/43073-a-guide-to-tags

Answer 1

[Jan]

Save this as file:

function c = func(t, y, s, q, w) 
c = [s * (y(2) - y(1)*(y(2)+1) - q*y(1)^2); ...   % α'(t)
     (-y(2) - y(1)*y(2) + y(3))/s; ...            % β'(t)
     w * (y(1)-y(3))];                            % γ'(t)

Then start the integrator like this:

s = 1;
q = 1; 
w = 0.1610;

y0 = [30 1 30];
t_span = [0 10];

[t, y] = ode45(@(t,x) func(t, x, s, q, w), t_span, y0);

See: http://www.mathworks.de/matlabcentral/answers/1971

[fxo]

Thank you for the help, worked very well!:)