Solving a differential equation numerically (in Octave, Matlab &c.)

[ejlflop]

Homework Statement

I have a second-order non-linear differential equation that I am trying to solve. So far I have decomposed it into a system of 2 first-order equations, and have (possibly) determined that it cannot be solved analytically. So I am trying to do a nice numerical approximation using GNU Octave (basically compatible with Matlab, so if you can do it with Matlab please can you help too :-D).

Octave needs the equation expressed as first-order DEs -- which I think I've done -- I'm just not sure how to go about doing the actual approximation.

Homework Equations

Original second-order DE:
$$\frac{d^2 \theta}{dt^2} + \frac{k}{m}\cdot\frac{d\theta}{dt} - g\cdot\sin\theta = 0$$

Note that $$k, m, g$$ are arbitrary constants (yes, $$g$$ is 9.81!)

The Attempt at a Solution

Substitute:
$$\frac{d\theta}{dt} = w$$

Hence:

$$\frac{d^2 \theta}{dt^2} = w\frac{dw}{d\theta}$$

$$w\frac{dw}{d\theta} + (\frac{k}{m}\cdot w) - (g\cdot\sin\theta) = 0$$

Attempt at Octave program to approximate it a bit:

function wdot = f (w, theta)
  g = 9.8
  k = 1
  m = 0.1
  wdot = (g*sin(theta) - (k/m)*w)/w
endfunction

theta = linspace(0, 20, 400);
y = lsode ("f", 1, theta);
plot (y, theta);

So this gives me a nice little graph, but obviously it's not a total solution -- that would involve computing the whole system of 2 DEs, which is what I don't know how to do!
Any help much appreciated, thanks.

 

[ejlflop]

No worries; I found a nice tutorial on how to do it in matlab, and adapted it for my own purposes.
If anyone has a similar problem, see this youtube video, it's very good:
https://www.youtube.com/watch?v=http://www.youtube.com/watch?v=fx3bl4oA_0U