INSTRUCTIONS FOR MATLAB CODES TO SOLVE ODES

The basic codes I've written to solve vector ODES are 
"ode.m", "dydt.m" and "f_RHS.m". You will need to 
download all of these files and place them in the 
same directory where you run matlab. 

The only code that you will need to open and change
yourself is "f_RHS.m". It is the program which 
contains the system of ODE's you want to solve. 
The sample program that is available here is for 
the pendulum equation. It looks like this:

function [dydt] = f_RHS( y, t ) 
dydt(1) = y(2); 
dydt(2) = -sin(y(1)); 

You can modify this file for any system of ODE's 
with any number of variables. Don't change the 
first line, but just write your equations into 
the following lines in the same format as above. 
[Don't forget the semicolons ";" at the end 
of each line! These prevent matlab from writing
the lines on the screen every time they're called
in the program.]

After entering your system of equations into the
file "f_RHS.m" you can solve the system by opening
matlab and running "ode.m". To do so, you must 
first enter the following information:

D, the dimension of the vector system 
y_0, the 1xD row vector of initial data 
t_0, the initial time
t_f, the final time
N_s, the number of integration steps
ord, the order of integrator (1=Euler,2=Heun,4=Runge-Kutta)
fig, the request for figures (1=yes, 0=no)

You can enter these at the matlab prompt. I'll 
illustrate this for the pendulum equation. Let's 
integrate it over the time interval from 0 to 5
with time-step h=0.01, requiring 500 steps, using
4-th-order Runge-Kutta. Let's take initial data 
y(0)=2.0, y'(0)=1.0. To do this you should type 
as follows at the matlab prompts:  

>> D=2;
>> y_0=[2.0,1.0];
>> t_0=0;
>> t_f=5;
>> N_s=500;
>> ord=4;
>> yes=1;

This enters all of the required data. Finally, type

>> [ y,t ] = ode(D,y_0,t_0,t_f,N_s,ord,yes) 

Try it! The code will produce a plot of the 1st component
y_1 vs. time t. If you hit 'enter', it will then produce 
a plot of the 2nd component y_2 vs. time t, and so on 
successively to the Dth component y_D vs. time t. Finally,
for D=2 (and only then), if you hit 'enter' a phase plot 
of the solution in the y_1y_2-plane will appear. The 
matlab prompt ">>" will then reappear. If you didn't wish 
any of these figures, then you could just have set fig=0 
at the beginning. In any case, the numerical solution 
will be stored when the program finishes and returns to 
the matlab prompt. The solution is written into a file 
called "data", where the first column is the time, the 
2nd column is y_1, the third column y_2, and so forth. 

If you want to solve the same ODE again with different
data, you can just enter the new data and run "ode.m" 
again. For example to solve the same problem as above
with h=0.001 or with 5000 time steps and with no figures, 
you would just type 

>> N_s=5000; fig=0;

All of the other data remains the same. Then type again

>> [ y,t ] = ode(D,y_0,t_0,t_f,N_s,ord,fig)

as before. Note: when you run "ode.m" again the file "data" 
will be overwritten with the new solution! If you wish 
to keep the data, you'll need to move it to another file 
before you run the code again. You can solve as many 
different cases as you like. If you wish to solve a 
new system of ODE's, then you must go into "f_RHS.m" 
and write in your new equations and run again. 

That's it!!!

By the way, matlab has its own ODE solvers which are 
much more sophisticated and accurate than the simple 
algorithms (Euler, Heun, 4th-order Runge-Kutta) that 
are implemented in my code. For example, one of these 
is "ode45". It is an adaptive Runge-Kutta-Fehlberg 
method, and a good "first try" for most problems. For 
information about it, just type 

>> help ode45

at the matlab prompt, or else go to this on-line URL:

 http://www.mathworks.com/access/helpdesk/help/techdoc/ref/ode15s.shtml

If you like, you can run it or other matlab ODE solvers
and compare results for speed and accuracy with the simple 
ones (Euler, Heun, 4th-order Runge-Kutta).