Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have to solve numerically this set of Ordinary Differential Equations

$$ \frac{dx_1}{ds} = \frac{1}{x_1} \left[x_2 \left(a + \frac{x_2}{s}\right)-\alpha x_1 z\right]$$

$$ \frac{dx_2}{ds} = -\frac{1}{x_1} \left[x_2 \left(b + x_2\right)z\right] - \left(c+\frac{x_2}{s} \right)$$

where $z(s) = \sqrt{x_1^2 +\left(b + x_2 \right)^2}$. The equations are time-independent, $s$ is a distance, $a(s)$, $b(s)$ and $c(s)$ are constant vectors, i.e. $f(s)$; and $\alpha$ is a constant. I've tried to solve using Matlab's ODE45 to no avail because perhaps the fact that the constants are vectors.

Here's the code and some vector examples:

global a b c d
y = ode45(@odeeqns, s, [1 2],[0 0]) % -> I'm not sure about these because I'm working with distance, not time.

function dyds = odeeqns(s,y)
   global a b c d
   dyds(:,1) = 1./y(1) .*(y(2) .* (a + y(2)./ s) - alpha .* y(1) .* sqrt(y(1).^2 + (b + y(2)).^2));
   dyds(:,2) = - (y(2) .* (b + y(2)) .* sqrt(y(1).^2 + (b + y(2)).^2))./y(1) - y(2)./s - c;
   plot(dyds);hold on

alpha = .0167;
a = [0.0000 0.2985 1.4973 2.4266 2.7838 2.7917 2.6397 2.4295 2.2094 2.0000];
b = [0.0000 0.0298 0.2246 0.4853 0.6960 0.8375 0.9239 0.9718 0.9942 1.0000];    
c = [0.0299 0.2615 0.9764 1.4490 1.5680 1.5098 1.3870 1.2499 1.1188 1.0058];    
s = [1 2 3 4 5 6 7 8 9 10]; % km

This is the error message:

??? Error using ==> odearguments at 113
YPRIME must return a column vector.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> windfield_powell at 33
[T,sols] = ode45(@yprime, [1 2],[1 1]);

Any help please? Thanks!

share|cite|improve this question
I don't really understand the "type" of your objects. When you say vectors you mean $a(s)$, $b(s)$, $c(s)$ are real numbers? And why do you say that the equations are time independent? ($s$ does appear on the right.) – Sebastien B Feb 11 '13 at 22:31
$s$ is a distance; yes, $a(s).$ are functions (expressed as vectors in Matlab) defined over distance $s$. – Oliver Amundsen Feb 11 '13 at 22:36
Ok, I think I understood your problem. Indeed I think it will be a problem to use ode45. What about coding yourself the explicit euler method or runge-kutta45? Then you can use your vectors. – Sebastien B Feb 12 '13 at 10:55
Do you maybe also have a refenrece to the exercise? Maybe a paper? – k1next Feb 12 '13 at 17:41
nope, did you see the update? – Oliver Amundsen Feb 12 '13 at 17:41

You can try something like that, first define a function

function dx=f(s,x)

a=1.2; b=2.3; c=3.4; alpha=4.5;





and then use ode45 with vectors:

[T,Y] = ode45(@f,[3 12],[1 1])


I hope it helps although perhaps you don't want to define the parameters $a$, $b$, $c$, $\alpha$ inside the function $f$

share|cite|improve this answer
The code is helpful, but remember that a, b and c are vectors of constants, not a single constant.... how should I implement that? – Oliver Amundsen Feb 11 '13 at 23:31
@OliverAmundsen can't you implement a function for a? as a = f(s)? – k1next Feb 12 '13 at 11:19
@macydanim, I think that the constant vectors $a$, $b$, and $c$ are already expressed as functions of $s$, in that they're vectors of the same dimensions. That being said, I also tried to integrate the ODE symbolically with Matlab syms, but they're to complicated... any other suggestions? – Oliver Amundsen Feb 12 '13 at 16:02
Just for my comprehension. The ODE system is the one given for $dx_1 /ds$ and $dx_2/ds$ as in the beginng of your question. But actually it isn't $a$ but $a(s)$ and the same with b,c. But these functions are explicitly know? So you could rewrite the ODE system in terms of $x_1,x_2,s$ without using $a(s),b(s),c(s)$ ? Could you give us the functions for $a,b,c$ ? Could you provide your Matlab code? That would help a lot! – k1next Feb 12 '13 at 16:28
Ah maybe I understood something. Do you have prior to running the ODE45 already the values for $a,b,c$ ? Assuming we want to evalute $x_1,x_2$ at $s=0:0.1:1$ then you already have a vector $a$ with $length(a)=10$? And you would like to hand that vector to the function that should be integrate and access the correct values, like for $s=0.2$ you want to access the third element of $a$ as $a(1)$ corresponds to $s=0$ and $a(2)$ to $s=0.1$ ,respectively? Or am I on a complete wrong course? – k1next Feb 12 '13 at 16:35

If you have a function then just add dyds=zeros(size(y))

Looks like this:

function dyds = odeeqns(s,y)


It should work.

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.