The following code doesn't seem to work - it keeps running with no result (RK4 method)

There are 3 differential equations as defined in 'f' below separated by semi-colon. The time span is from 0 to 84 by 7.

f = @(t,x) [388-0.009*x(1)-0.0000000361*x(1)*x(3); 0.0000000361*x(1)*x(3) - 0.18*x(2); 50000*x(3) - 23*x(3)];
tspan = 0:7:84;
[t,xa] = ode45(f,tspan,[43100 0 0.0033]);

[Walter Roberson]

Could you post the original equations?

[KC]

x' = λ - ρx - βxz;
y' = βxz - δy;
z' = py - cz;
x0=43100; y0 = 0, z0 = 0.0033, λ = 388, ρ = 0.009 δ = 0.18, p = 50000, c = 23, β=3.61e-8

[Walter Roberson]

It appears that you need to add β = 0.361e-7 to match the f you posted.

Answer 1

[Marc]

You figure out your mistake yet??

Your code is just "fine"... But your function as you have it above is "wrong".

Take a look at your third equation. Change that first x(3) to a x(2) and you will be in business.

No issues here using ODE45 or ODE15s. Also, you don't need to set tspan as such... [0:7:84] You can simply use [0 84] and you get a smoother answer. This took well under a second for both 45 and 15s....

Interestingly, you should try both and compare. ODE45 gives you a solution in ~0.1-0.2 seconds (depending on your machine) and the solution is ~2400 points. 15s gives you an answer in under 0.1s and the solution is 163 points. Just goes to show you how much more efficient ODE15s can be with stiff problems.

[KC]

Thanks Marc!

Answer 2

[John D'Errico]

This is surely what would be called a stiff problem. It will cause ODE45 (or any non-stiff solver like it) to run very poorly. The first clue to the problem being a stiff one is the fact that ODE45 DID work poorly. A better clue is the wide disparity of magnitudes of the coefficients in that ODE. That tends to create stiff systems.

So instead, try one of the stiff solvers. All of the stiff solvers have names that end in s, thus ode23s, ode15s.

[KC]

I just tried using both ode23s and ode15s, but I get many warning messages like these: > In ode23s (line 379) Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.484584e-246. Warning: Failure at t=7.014056e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.491893e-14) at time t. I am wondering if my code is correct?