How can I solve ordinary differential equations in MATLAB?

Matlab can numerically solve Ordinary Differential equations using 2 methods.

ODE23 uses 2nd and 3rd order Runge-Kutta formulas
ODE45 uses 4th and 5th order Runge-Kutta formulas

What you first need to do is to break your ODE into a system of 1st order equations. For instance,

        ..    .    2
        X  + AX + K X = 0       <= Damped harmonic oscillator

can be written as,

        .
        Y1 = Y2
        .              2
        Y2 = - A Y2 - K Y1
                 .
(Let Y1=X and Y2=X)

Now, you need to write a matlab function that takes Y1, Y2, and time as arguments and returns Ydot1 and Ydot2. To do this, create a file using emacs or your own favorite editor that contains the following:

    function [Ydot] = myode(t,Y)

    % Note: Y(1) => Y1 and Y(2) => Y2
    % t is for time.. if your equations do not contain t, you still have
    % to put a "t" in the function declaration above.

    A = 1.3333;    % Whatever values you want to choose
    K= 2.132;

    Ydot(1) = Y(2);
    Ydot(2) = -A*Y(2)-K^2*Y(1);

    Ydot = Ydot';  % This makes Ydot into a column vector

Call this file "myode.m" and save it in a place where matlab knows to look for it (either the directory from which you started matlab or your own /mit/username/matlab/ directory).

Now, the moment we have been waiting for: let's solve this ODE. At the matlab prompt, type:

	 >> [T,Y]=ode45('myode',[0 10],[0 1]);

Note:

        'myode': Name of the function you wrote above
        [0 10] : Integrate from t=0 to t=10 (seconds)
        [0 1]  : Initial conditions for your system.. [1 0]
                 corresponds to system where,

                        X(t=0) = 1
                        .
                        X(t=0) = 0

                 Which you can think of as a pendulum that has
                 been pulled back and is being let go at t=0.
                 (in this case the pendulum is damped by the A term)

To see the results, type:

        >> plot(T,Y)

or,

        >> plot(T,Y(:,1),'-',T,Y(:,2),'--');

Here, the solid line corresponds to X (displacement) and the dashed line to Xdot (speed).

Your coefficients to your equations do not have to be constant. They can depend on time. For instance if your oscillator is damped and driven:

        ..    .    2
        X  + AX + K X = B cos(wt)

Your m-file will look something like this:

    function [Ydot] = myode2(t,Y)

    % Note: Y(1) => Y1 and Y(2) => Y2
    % t is for time.. if your equations do not contain t, you still have
    % to put a "t" in the function declaration above.

    A = 1.3333;    % Whatever values you want to choose
    K= 2.132;
    B=.2;
    w=2;

    Ydot(1) = Y(2);
    Ydot(2) = -A*Y(2)-K^2*Y(1)+B*cos(w*t);

    Ydot = Ydot';  % This makes Ydot into a column vector

You might try experimenting with the coefficients to see what happens.

    >> [T,Y]=ode45('myode2',[0 10], [0 0]); plot(T,Y);